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ABSTRACT 

A simple ID dynamical model of thermally driven disc winds is proposed, based on 
the results of recent, 2.5D axi-symmetric simulations. Our formulation of the disc 
wind problem is in the spirit of the original Parker (1958) and Bondi (1952) problems, 
namely we assume an elementary flow configuration consisting of an outflow following 
pre-defined trajectories in the presence of a central gravitating point mass. Viscosity 
and heat conduction are neglected. We consider two different streamline geometries, 
both comprised of straight lines in the {x, z)-plane: (i) streamlines that converge to 
a geometric point located at (x, z) = (0, —d) and (ii) streamlines that emerge at a 
constant inclination angle from the disc midplane (the x-axis, as we consider geo- 
metrically thin accretion discs). The former geometry is commonly used in kinematic 
models to compute synthetic spectra, while the latter, which exhibits self-similarity, 
is likely unused for this purpose, although it easily can be with existing kinematic 
models. We make the case that it should be, i.e. that geometry (ii) leads to transonic 
wind solutions with substantially different properties owing to its lack of streamline 
divergence. Both geometries can be used to complement recent efforts to estimate pho- 
toevaporative mass loss rates from protoplanetary discs. Pertinent to understanding 
our disc wind results, which are also applicable to X-ray binaries and active galactic 
nuclei, is a focused discussion on lesser known properties of classic Parker wind so- 
lutions. We find that the parameter space corresponding to decelerating Parker wind 
solutions is made larger due to rotation and leads instead to disc wind solutions that 
always accelerate after the bulk velocity is slowed to a minimum value. Surprisingly, 
Keplerian rotation may allow for two different transonic wind solutions for the same 
physical conditions. 

Key words: accretion discs - hydrodynamics ~ planets and satellites: atmospheres 
- protoplanetary discs - stars: winds, outflows 



1 INTRODUCTION 

The classic Parker model has served as a paradigm wind so- 
lution for over half a century now. First developed for the 
Sun as a model of the solar wind, it shows the essential 
features of one-dimensional (ID), steady state wind models, 
namely that transonic solutions typically involve a transition 
through a critical point and have an X-type solution topol- 
ogy. The analytic solutions to both the original isothermal 
(with adiabatic index 7 = 1) Parker wind model (Parker, 
1958) and its polytropic fluid (1 < 7 < 5/3) extension 
(Parker, 1960) serve a dual purpose. On the one hand, they 
are beneficial for obtaining insight into the theory of outfiows 
in general, as well as for gaining intuition into the subtleties 
that arise when solving wind equations analytically (for an 
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in depth perspective see Konigl & Salmeron 2011 and ref- 
erences therein). On the other hand, Parker wind solutions 
have proven useful for assessing the accuracy of numerical 
simulations (e.g., Keppens & Goedbloed 1999; Font et al. 
2004; Tian et al. 2005; Stone & Proga 2009). To that end, 
one goal of this paper is to present a simplified dimension- 
less formulation of Parker winds and to provide formulae 
commonly used for numerical testing purposes. At the same 
time, we address certain aspects of the polytropic Parker 
problem that have been a source of confusion in the liter- 
ature. Specifically, we clarify the properties of spherically 
symmetric Parker winds in the range 3/2 < 7 < 5/3 and 
the corresponding range of 7 when angular momentum is 
added to the problem. 

The primary focus of this paper is to present solutions 
to Parker-like winds emanating as a biconical flow, the geom- 
etry commonly used to model accretion disc winds. Parker 
winds have been instrumental in uncovering other physical 
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processes that can drive winds in stars, and led to the de- 
velopment of both line-driven (Castor et al. 1975, hereafter 
CAK) and magneto-centrifugally driven (beginning with the 
solution of Wcbcr &: Davis 1967) wind theory. The current 
state of the art in stellar wind theory owes much of its de- 
velopment to the systematic assessment of how the inclusion 
of various physical terms and geometrical effects in the hy- 
drodynamic equations alters the solutions of Parker winds. 
Studies of disc winds stand to benefit from rigorously re- 
peating this procedure using counterpart, ID analytical disc 
wind models. 

Developing concrete baseline models analogous to 
Parker winds has proven to be a difficult task. A major road- 
block has been the uncertainty in the streamline geometry, 
i.e. the actual trajectory traversed by gas flowing out from 
the disc, as well eis in the gravitational potential along these 
streamlines. Another obvious and related difficulty is posed 
by the fact that accretion discs span many more orders of 
magnitude in physical size than do stars, and they can host 
radically different, spatially and temporally variable, ther- 
modynamic environments. Indeed, the outer radius of an 
accretion disc ranges from parsec scales for active galactic 
nuclei (AGN) down to within 1 AU for some circumstellar 
discs and the diverse physical conditions permit anything 
from infrequent outbursts to highly relativistic, steady jets. 
It should come as no surprise then, that despite clear obser- 
vational evidence of outflows from many systems, identifying 
the actual driving mechanisms, as well as determining the 
wind geometry, remains a challenge. 

Studies of disc winds therefore rely heavily on kinematic 
models in order to quickly explore the parameter space with- 
out assuming a particular driving mechanism. For exam- 
ple, kinematic models have been employed to produce syn- 
thetic spectra for cataclysmic variables (CVs), systems in 
which even key properties such as the geometry, ionization 
structure, and mass-loss rates remain difficult to constrain 
(e.g., Noebauer et al. 2010 and references therein). Early 
kinematic models assumed spherically symmetric outflows 
for simplicity (Drew & Verbunt 1985; Mauche & Raymond 
1987). The consensus picture of a biconical mass outflow 
originating from the inner disc was born out of the observed 
characteristics of resonance lines in CVs (Cordova & Mar 
son 1985, Drew 1987). This geometry was developed into 
a robust kinematic model by Shlosman & Vitello (1993), 
who calculated the ionization structure of CV disc winds 
and solved a radiative transfer problem in lines using the 
Sobolev approximation. Their kinematic model allowed for 
an arbitrary amount of streamline divergence. 

Knigge et al. (1995) developed a different (Monte Carlo) 
code to solve the radiative transfer exactly. Their choice of 
wind geometry is one instance of what we refer to as the Con- 
verging model, in that the streamline divergence is just such 
that all streamlines converge to a geometric point located a 
distance d below the disc, as illustrated in Figure lb. The 
Converging model, which has been called the "displaced- 
dipole" model by others, has been used in conjunction with 
sophisticated radiative transfer simulations to model accre- 
tion disc spectra from massive young stellar objects (Sim et 
al. 2005), active galactic nuclei (Sim et al. 2008), CV disc 
winds (Noebauer et. al 2010), classical T Tauri stars (Kuro- 
sawa et al. 2011), and young intermediate-mass Herbig Ae 
stars (Grinin & Tambovtseva 2011). Typically, these sim- 




c. CIA Model 




Figure 1. Pictorial representation of the streamline geometry 
addressed in this paper. Neighboring streamlines diverge from 
each other in the (a) Parker and (b) Converging models, whereas 
in (c), the Constant IncHnation Angle (CIA) model, there is no 
adjacent streamline divergence. 

ulations use Monte Carlo procedures that can account for 
nearly all of the prominent resonance lines and thereby ac- 
curately calculate the ionization balance of the wind. The 
Converging model has even been employed to calculate the 
neutron structure of neutrino-heated MHD disc winds asso- 
ciated with gamma-ray bursts (Metzger et al. 2008). 

In this paper we develop a simple dynamical disc wind 
model that amounts to a generalization of the Parker model. 
Rather than positing a velocity law as is done for kinematic 
models, the purpose of a dynamical model is to impose the 
physical conditions and solve for the wind velocity as a func- 
tion of distance along a streamline. This necessarily requires 
identifying a driving mechanism, i.e. a heating source in the 
case of thermally driven winds. Much of the groundwork 
theory for the source of heating was laid down by Begelman 
et al. (1983, hereafter BMS83), who showed that Compton- 
heated coronae are qualitatively the same for both quasars 
and X-ray binaries. That is, both galactic X-ray sources and 
the inner regions of AGN are expected to be heated via ir- 
radiation from a central X-ray source to high enough tem- 
peratures that thermal expansion alone gives rise to a disc 
wind. 

The hydrodynamic formulation of BMS83 established 
that it suffices to estimate 2D global wind properties with a 
ID model that captures the essential physics. Indeed, many 
predictions given by BMS83 were later confirmed by fol- 
lowup works that focused on the inherently two-dimensional 
radiative transfer problem (e.g., Ostrikcr et al. 1991, Woods 
et al. 1996, Proga & Kallman 2002) . Of special interest here 
is the work by Woods et al. (1996), who added to the ba- 
sic theory of BMS83 based on the outcome of their time- 
dependent, 2.5D simulations of thermally driven winds from 
AGN heated by Compton as well eis non-Compton processes 
such as photoionization and line-cooling. They provided an 
improved formula for the mass flux density as a function 
of disc radius and presented a detailed study of the flow 
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2 PARKER-LIKE DISC WINDS 
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Figure 2. Flow streamlines that resulted from the time- 
dependent, hydrodynamical simulation of a thermally driven wind 
(Luketic et al. 2010). The 2-axis is the rotation axis, while the x- 
axis is the disc midplaue. Streamlines at x > ^Ric i where Rjc is 
the Compton radius, are self-similar. This figure gave motivation 
for the CIA model. 



topology and sonic surfaces for various spectral energy dis- 
tributions. 

Both the results of Woods et al. (1996) and those of the 
more recent 2.5D time-dependent simulations of a thermally 
driven wind carried out by Luketic et al. (2010, see Figure 2 
here) indicate that the streamline geometry is rather simple, 
displaying two distinct flow patterns. Moreover, their results 
suggest that the Converging model may not be well-suited 
for sampling the entire wind, but rather only the inner por- 
tions of it. The outer portion is better approximated by a 
model in which streamlines emerge at a constant inclination 
angle to the midplane (hence the name, the CIA Model - 
see Figure Ic).^ It is our intention to study how this dif- 
ference in geometry affects the hydrodynamics independent 
of the explicit heating mechanism taking place; we merely 
assume that the boundary of the flow (the disc midplane) 
has been heated to a high enough temperature to drive a 
thermal wind. 

This paper is organized as follows. Background theory 
pertaining to the hydrodynamic foundations of wind theory, 
sources of thermal driving, and our choice of streamline ge- 
ometry is presented in §2. Our formulation is given in §3. We 
make a comparison with simulations and present our results 
in §4. We discuss subtle aspects of the classic Parker prob- 
lem that arc relevant to disc winds and reveal the effects of 
adding rotation in §5. Finally, we summarize our findings 
and open questions in §6. 

^ The kinematic model used by Shlosman & Vitello (1993) can 
accommodate CIA streamlines by setting 6min = dmax- 



Aside from stellar winds, spherical Parker winds have been 
utilized to model AGN winds (e.g., Everett & Murray 2007) 
and protoplanetary disc winds (e.g., Adams et al. 2004, 
Gorti & HoUenbach 2009, Owen et al. 2012). They are also 
very useful for conceptualizing the atmospheric escape pro- 
cess from planets (e.g.. Tucker et al. 2012 and references 
therein). By relaxing the assumption of spherical symme- 
try, we call the outflow traversing an azimuthally symmet- 
ric streamline geometry a Parker-like disc wind. The lat- 
ter retains the same notion of a clcissic Parker wind, which 
we review in §2.1 before discussing concepts specific to disc 
winds. 



2.1 Static Atmospheres vs. Parker Winds 

Parker winds describe a highly idealized fluid phenomenon: 
the steady state, spherically symmetric, hydrodynamic es- 
cape of an ideal gas with zero viscosity. With these simpli- 
fications, and further supposing a fiuid composed of only 
one species of gas, the Eulerian equations become analyt- 
ically tractable. A complete isothermal Parker wind solu- 
tion consists of the density p and velocity as a function 
of radius r. Polytropic Parker wind solutions also include a 
temperature profile T{r). The family of isothermal (7 = 1) 
transonic Parker wind solutions depends solely on one pa- 
rameter, which is often called the hydrodynamic energy pa- 
rameter (HEP) and is defined as 



(2.1) 



■ykToVo 

Here, G is the gravitational constant, M, is the central 
body's mass, nip the proton mass, fi is the mean molecular 
weight, k is Boltzmann's constant, and To is the tempera- 
ture at the base of a streamline, i.e. at the boundary radius 
To. A polytropic equation of state (EoS) introduces a second 
dependence on the adiabatic index 7. (The choice to incor- 
porate 7 into the HEP was not made in the early papers 
on the solar wind, but it makes the dimensionless equations 
less cluttered.) 

Physically, Parker wind solutions model atmospheric 
coronae in hydrodynamic equilibrium just as a familiar baro- 
metric law, which yields the variation of density with radius, 
models an atmosphere in hydrostatic equilibrium. Indeed, 
an isothermal barometric law (albeit one that accounts for a 
varying gravitational potential) can be considered the trivial 
Parker wind solution with v{r) = 0, T{r) = To = constant, 
and 



P{r) = Poexp [Xo{ro/r) - Xo] ■ 



(2.2) 



where po is taken to be the known density at a reference 
radius Vo- Equation (2.2) is most commonly recognized as 
the solution to the equation of hydrostatic equilibrium, 

dP _ GM,p 
dr 



(2.3) 



for the special case of an ideal gas with a pressure P = 
pkTo/nipfi. However, equation (2.2) is readily seen to be 
the density profile found by taking the limit of a Parker 
wind solution with a slowly expanding atmosphere {v — ^ 0) 
and a small mass-loss rate (corresponding to an everywhere 
subsonic solution - see §3.7 for more discussion). Recovering 



© 0000 RAS, MNRAS 000, 000-000 



4 Waters & Proga 



a static atmosphere from a slowly expanding one hints at a 
correspondence between thermal escape processes in kinetic 
theory and fluid dynamics. 



2.1.1 Hydrodynamic Escape 

Parker winds capture the simplest example of a more gen- 
eral thermal escape process characterized by hydrodynamic 
escape. Modern, more realistic models can account for addi- 
tional physical processes and non-spherical geometries, but 
the underlying hydrodynamic, thermal escape mechanism is 
effectively isolated by Parker winds. 

The thermal escape process from a static atmosphere 
(i.e. evaporation or Jeans escape) is governed by kinetic 
theory and sometimes referred to as hydrostatic escape to 
distinguish (and emphasize) its relation to hydrodynamic 
escape (e.g., Seager 2010, pg. 448). A parameter common to 
both the kinetic theory and fluid dynamics approaches to de- 
riving a barometric law for a static (i.e. slowly evaporating) 
atmosphere is what we call the thermal energy parameter 
(TEP), 



l$l 



T{r)/To' 



(2.4) 



where $ = —GM,/r is the gravitational potential energy 
per unit mass and Cs = ^/ 'ykT / jimp is the adiabatic speed 
of sound. The TEP is by definition a measure of the thermal 
energy of gas at every location in a central gravitational field, 
and the second equality permits us to interpret the HEP as 
just the TEP evaluated at some reference level ro, at which 
the temperature is To. The magnitude of the HEP at this 
level, as well as the asymptotic value of the TEP, governs 
which approach, fluid or kinetic, better models the escape 
process (for concrete examples, see Tucker et al. 2012; see 
also Owen &: Jackson 2012). 

Transonic hydrodynamic escape is associated with the 
important property that the pressure tends to zero asymp- 
totically (Parker 1958, 1960, 1965). That static atmospheres 
can lack this property in essence provided the physical ba- 
sis for Parker's original transonic solar wind solution, as 
the hydrostatic conduction model of Chapman (1957), the 
model of the extended solar corona that Parker's model su- 
perseded, implicitly featured a non-vanishing pressure at in- 
finity.^ Parker reasoned that unless the pressure vanishes at 
infinity, a static atmosphere can only truly be held static if 
there is a finite inward pressure exerted on it at large radii. 
In the case of the Sun, Parker pointed out that the vacuum- 
like conditions of the interstellar medium cannot possibly 
provide the necessary back pressure to keep the Sun's atmo- 
sphere in hydrostatic equilibrium (although see Velli 2001 
and references therein for the extent to which this argument 
holds). 

More rigorously, we can exploit the TEP to identify a 



^ Indeed, the isothermal barometric law of equation (2.2) has a 
non- vanishing pressure at infinity as the density in equation (2.2) 
at r = GO is poexp(— Ao), and for barotropic flow, the pressure is 
a function only of density. See Chamberlain (1963) for a detailed 
discussion of this breakdown of a barometric law using a kinetic 
theory approach. 



threshold temperature decline that determines whether or 
not an atmosphere can be held in hydrostatic equilibrium. 
Integration of equation equation (2.3) over a non-isothermal 

atmosphere that extends from To (where p = po and T = To) 
to some radius r can be written (Parker 1965), 

r(r') 



p(r)T(r) = poToexp 



- dr 



(2.5) 



where wc have taken the pressure as P(r) = Cs{r)^ p{r). 
In order for equation (2.5) to describe a static atmosphere 
surrounded by vacuum, the density must vanish at infinity, 
implying that the integral inside the exponent must be di- 
vergent at large r. Conversely, a Parker wind is the steady 
equilibrium state of an atmosphere if the integral is conver- 
gent. In terms of the TEP, p vanishes at infinity if t is an 
increasing function of r, while the density tends to a finite 
value if r(r) is decreasing. Physically, therefore, a spherical, 
isolated static atmosphere is possible only if the magnitude 
of the gravitational potential energy of the gas outweighs its 
thermal energy at large radii. For the critical case in which 
these energies are in balance, i.e. when r(r) is constant, we 
see from equation (2.4) that the temperature profile satisfies 
T{r)/To = ro/r (provided Ao does not vary with r) and from 
equation (2.5) that the integral diverges logarithmically. 

It seems to have been overlooked previously that there 
are no transonic Parker wind solutions with a l/r tem- 
perature dependence, which occurs when 7 = 3/2 (sec 
Appendix D). Importantly, 7 = 3/2 is the critical adiabatic 
index that divides the behavior and solution space of tran- 
sonic polytropic Parker winds. For 7 < 3/2, transonic Parker 
winds are accelerating, while for 7 > 3/2, they are decel- 
erating. As summarized in Tabic 1, the decelerating wind 
regime permits isolated hydrostatic solutions (which have 
dr/dr > 0), while only Parker winds can have a vanish- 
ing pressure at infinity for 1 ?J 7 < 3/2. It is clear from 
Table 1 that the parameter space, (Ao,7), leading to spher- 
ically symmetric transonic Parker wind solutions is coupled 
in a simple way. In §5, these HEP bounds are generalized to 
account for the effects of rotation. 



2.2 The Applicability of Parker-like Disk Winds 

Magncto-ccntrifugally driven winds are often invoked as can- 
didate mechanisms for explaining outflows from accretion 
discs. In systems or regions of systems where magnetic forces 
might be dynamically unimportant, thermal driving is a 
likely contender (e.g., Proga 2007 & references therein). Just 
as Parker winds can be useful for modeling outflows from any 
spherical astrophysical body thought to be hot enough to ex- 
hibit a non-explosive, thermal expansion of gas, the Parker- 
like disc winds addressed in this paper can be used to model 
thermally driven winds from the coronae of accretion discs 
associated with AGN, X-ray binaries, unmagnetized proto- 
stellar discs, and protoplanotary discs. Due to the diversity 
of physical scales in these systems, a preliminary step for 
constructing a ID model is to identify a characteristic rar 
dius for invoking thermal driving, in order to calculate the 
HEP. First, it is worth emphasizing that the escape velocity 
for discs varies with distance To along the disc midplane, so 
the HEP must be considered a function of ro. In this regard 
there is an intrinsic difference between ID disc wind models 
and ID spherical wind models. Namely, for a given mass M* 
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Table 1. Parameter Space of Spherically Symmetric Parker Wind Solutions. 



Polytropic 


Permitted 


TEP 


Hydrostatic^ 


Transonic 


Index 


HEP Range 


Behavior 


solutions? 


solutions? 


7 = 1 


[2,oo] 


dr/dr < 


No 


Yes 


1< 7 < 3/2 


[2,1/(7-1)] 


dr/dr < 


No 


Yes 


7 = 3/2 


Ao = 2 


dr/dr = 


Yes 


No 


3/2 < 7 < 5/3 


[1/(7-1), 2] 


dr/dr > 


Yes 


Yes 


7 = 5/3 


[1.5,2] 


dr/dr > 


Yes 


No 



1'This refers to 'isolated' hydrostatic solutions, i.e. those with a vanishing density at infinity. 



and characteristic launching radius To for a star or planet, 
varying the HEP samples different temperatures of the stel- 
lar corona or planetary exosphere. Meanwhile, varying the 
HEP for a given central object mass for disc winds corre- 
sponds to altering either the temperature at a fixed distance 
along the midplane or the distance at a fixed temperature - 
or both. 



2.2.1 The Thin-Disc Assumption 

The scale height of an isothermal corona is given by H = 
-y/2/Ac, ro, implying that our models can only be applied to 
regions of a flared disc where Ao >> 2, as we are implic- 
itly supposing that H/ro « 1. In order for our disc wind 
solutions to not be restricted to Ao >> 2, the sound speed 
within the disc must be considered separate from the sound 
speed at the base of the wind. In other words, we imagine 
a cold, thin disc that acts as a reservoir of material capable 
of sustaining a wind. Models of the internal disc structure 
show it to be complex and turbulent (e.g., Balbus & Hawlcy 
1998; Miller & Stone 2000; Proga & Bcgclman 2003; Turner 
et al. 2003; Hirose et al. 2006; Blaes et al. 2007; Krohk et 
al. 2007). Since there is no analytic model for this dynamic 
internal structure of the disc, wc cannot incorporate this 
complexity into our treatment. Hence, the gas in the disc 
need not 'match' onto the base of the streamlines of the 
heated surface gas that forms the outflow, as the situation 
is analogous to the photosphere-corona transition. 



2.3 Sources of Heating for Thermal Driving 

The two simplest heating mechanisms believed to be capable 
of launching thermal winds in accretion discs are Compton 
heating and photoionization. Here we discuss how to approx- 
imate the different wind regimes identified by BMS83 using 
the two input parameters (Ao,7) of Parker winds. 



2.3.1 Compton Heating 

The relevant length scale for AGN and X-ray binary disc 
winds is the Compton radius, the radius where the gravita- 
tional and thermal pressures are equal: 

GM« _ GM^mpfx 



Ric = 



kTic 



(2.6) 



Here, cic is the isothermal sound speed for gas heated to the 
inverse Compton temperature Tic (which can be ~ 10® K 



depending on the spectrum of radiation) , defined by 

kTic = \ {hv) , (2.7) 



where IJiv) is the average photon energy from an isotropic 
radiation source of luminosity C, namely 



1 



hvCu dv. 



(2.i 



As discussed by BMS83, regardless of magnitude of the lu- 
minosity, at radii beyond Ric the gas cannot remain quasi- 
static; the corona is itself unbound and better described as 
a vigorous wind region. Woods et al. (1996) found the di- 
viding cutoff for weak outflows to lie at the smaller radius 
To ~ 0.1 -Ric (see also Proga & Kallman 2002). In terms of 
C = ro/Ric, the HEP is 



A<, = — 



1 



Tic 



(2.9) 



Depending on the luminosity, the wind regions to either 
side of ~ 0.1-R/c are further divided; BMS83 identified five 
solution regimes in all (see also Woods et al. 1996). Each 
has an associated mass flux density, determined by ^ and 
C/ Ccr, where Ccr is a critical luminosity defined by 



1 



(- 



8/i \cic J \'m 



(2.10) 



Here, £b is the Eddington luminosity, c the speed of light, 
and me the electron mass; Lcrj Ce < 0.1 for Tic ^ lO'^ 
K, so that thermal pressure dominates radiation pressure. 
Parker-like disc winds are applicable in the regions affected 
by gravity, which includes the two weak wind regions with 
^ < 1, labelled D and E by BMS83 with £/£cr < 1 and 
L/ Lor > 1, respectively, as well as the portion of the 'grav- 
ity inhibited' strong wind region, labelled C, with ^ > 1 and 
C/Ccr « 1. The remaining two regions A & B have high 
enough luminosities that gravity is dynamically unimpor- 
tant and adiabatic losses insignificant in the subsonic flow 
regime. 

Utilizing Parker-like disc winds in the context of Comp- 
ton heating amounts to a significant simplification of the 
theory developed by BMS83. However, our models may pro- 
vide an adequate approximation of the disc wind dynamics 
because the functional form for how the mass flux density 
scales with radius fo along the disc plane is identical with 
that found by BMS83 for both isothermal (7 = 1) and isen- 
tropic (7 = 5/3) flow. We can account for this agreement by 
contrasting the two approaches used to treat the thermody- 
namics. 
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The simplicity of invoking Parker winds resides in the 
use of a polytropic EoS {P oa p~' , where P is the pressure and 
p is the density), the conventional means for bypassing the 
heat equation when the source of heating is very complicated 
or poorly understood (see, e.g., Tsinganos & Trussoni 1990; 
Sauty et al. 1999; Mehani et al. 2004) - not the case with 
Compton-heated Corona, in which the thermodynamics can 
be conveniently handled via an entropy equation (BMS83). 
If it is assumed that no heat is transferred via conduction or 
viscous dissipation to or from outflowing gas, conservation of 
energy dictates that the entropy production is proportional 
to the heating rate, T. For optically thin gas heated to tem- 
peratures T > 10® the net heating and cooling rate is 
proportional to the difference AT = T — Tic (e.g., Krolik et 
al. 1981). We therefore see that when the heating rate is high 
throughout the entire subsonic wind region, so that there is 
a near balance of heating and cooling (AT = 0), then an 
isothermal Parker wind with T = Tic will be a good ap- 
proximation to BMS83's strong gravity, nearly isothermal 
region E (see §4.1.2 for details). 

In the opposite case of adiabatic, iscntropic flow (7 = 
5/3), the entropy production is zero. In the framework of 
BMS83, this can effectively occur when the heat-transport 
can altogether be ignored (F as 0), meaning that the heat- 
ing time-scales are long compared to the flow time-scales. 
More speciflcally, 7 = 5/3 applies to gas with no internal 
degrees of freedom that is heated to a high temperature 
T < Tic in, say, a thin layer above the optically thick disc, 
that from there expands outward, loses additional pressure 
support upon being slowed by gravity, and thereby adiabat- 
ically cools. In §4.1.2, we explicitly show that the functional 
form of the mass flux density for a 7 = 5/3 Parker wind is 
identical to the prediction given by BMS83 for their solution 
regime C. 



2.3.2 Photoionization Heating 

As discussed by BMS83, very similar physics underlies 
Compton and photoionization heating, albeit the cooling 
mechanism for the latter is significantly more complicated 
(line-cooling and recombination vs. inverse Compton). Due 
to these complications, recent analytical studies have in- 
voked a combination of numerical simulations and spheri- 
cal Parker wind solutions in order to estimate global mass 
loss rates from protoplanetary discs (e.g., Gorti & HoUen- 
bach 2009, Owen et al. 2012). With the caveat mentioned in 
§2.2.1, our models make it possible to move beyond a spher- 
ical wind boundary and analytically investigate Parker-like 
winds from the surface of the disc. The starting point for our 
models is to make an explicit comparison between equation 
(2.9) and the HEP for protoplanetary discs, 

Ao = Vg/ro, (2.11) 

where rg = GMt/cl is 'the gravitational radius', the dis- 
tance where the gas becomes unbound because the escape 
velocity from the disc is equal to the thermal velocity of 
the gas. For a constant temperature on the disc midplanc 
(appropriate for a disc surface heated by EUV radiation to 
~ W^K or Compton heated to To{^) ~ 10^ K), rg is con- 
stant, and we see that photoevaporative winds are qualita- 
tively similar to Compton heated winds in the sense that rg 



plays the role of Ric. In either case, Ao decreases as ^ due 
to the reduced escape speed. 

2.4 Disc Streamline Geometry 

There are two routes to take in regards to specifying a geom- 
etry when finding solutions to wind equations (e.g., BMS83; 
Tsinganos & Sauty 1992): (i) assign some trajectory to the 
flow emanating from the disc or (ii) self-consistently solve 
for that trajectory. The first entails that an expression be 
provided for the flow tube area A that enters the steady 
state continuity equation pvA = constant. By adopting the 
geometry of the models discussed in the introdxiction, we 
necessarily take route (i). As shown in Figure 3, our geome- 
try is comprised of streamlines that are straight lines in the 
(a;, «)-plane. By rotating these straight streamlines around 
the 2:-ax;is, the actual trajectory traversed by the gas as it 
rises above the disc can be visualized; it spirals about a cone 
that widens according to the inclination angle i. Besides be- 
ing obscrvationally motivated, a biconical flow area is the 
simplest possible choice, for the distance along a streamline 
I = l{x, z) can be used as the sole variable instead of seek- 
ing some relationship bctwoon the cylindrical coordinates 
X and z. Our coordinates are related by x = To -I- Icosi, 
z = Isin i, and r = \^r'i + P + 2lro cos i. CIA streamlines 
are self-similar, while Converging streamlines are not. We 
find A{1) for each configuration in §3.1. 

BMS83 and Fukuc (1989) also took route (i) by as- 
suming a flow conflguration. Fukue (1989) adopted an area 
function similar to that used by BMS83, but he did not 
self-consistently implement the polytropic EoS when he fol- 
lowed BMS83 in requiring that the wind be launched from 
rest from the disc midplanc. (No restriction is imposed on Vo 
for the entropy equation used by BMS83.) In our notation, 
BMS83 chose a generic area function aimed at parametrizing 
the streamline divergence: A{1) = [l + l/roY . The parameter 
q is constrained to lie between (vertical flow) and 2 (spher- 
ical flow). In §3.1, we show that the Parker and Converging 
models have q = 2, while the CIA model has q = 1. 

Alternatively, route (ii) can be followed, in which an 
attempt can be made to calculate A{1) as part of the solu- 
tion. This involves either solving the fully 2D, 2.5D, or 3D 
problem using numerical techniques or it requires making an 
extra assumption, such as self-similarity or force balancing. 
The latter route was taken by Takahara et al. (1989), who 
arrived at an expression for A{1) by assuming that the cen- 
trifugal force balances the component of the gravitational 
force jjcrpendicular to the flow velocity at every distance I 
along the streamline.^ As discussed by BMS83, this is a valid 
approximation regardless of the streamline trajectory close 
do the disc midplanc if angular momentum is conserved. It 
was soon pointed out by Fukue & Okada (1990) that Taka- 
hara et al. (1989) misreiircscnted the gravitational force in 
calculating this force balance.'' The correct expression for 

^ In the MHD literature, balancing forces perpendicular to the 
streamlines leads to the Grad-Shafranov (or transfield) equation. 
** Taltahara et al. (1989) used the x-component of the gravi- 
tational force rather than the component perpendicular to the 
streamline, to arrive at a self-similar streamline trajectory given 

hy z = x-\J (x/to)'^/^ — 1. 
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Figure 3. The biconical outflow of our 'global' disc wind model is shown schematically on the left. On the right, we depict our coordinate 
system and the streamline geometry of the CIA and Converging models. In the fourth quadrant, we illustrate the geometry used to arrive 
at an expression for A(l), which can be visualized as the area swept out around the z-axis by any two neighboring streamlines at a fixed 
distance I. This notion becomes ex3x:t in the diflferential limit of closely spaced streamlines. 



A{x,z), obtained by Fukue & Okada (1990), is significantly 
more complicated; it suffices to consider the shape of their 
streamlines in the {x, 2)-plane. At every footprint distance 
To, the local streamline is found from 



To + cz 



1 + 



2z 



To + CZ 



(2.12) 



where r — \/ x"^ ^ z^ is the spherical position coordinate. The 
parameter c is an integration constant that resulted from the 
correct treatment of the force balance, which entailed solv- 
ing an ordinary differential equation for dx/dz. Hence, c is 
related to the slope - the opening angle of the streamline 
- and defines a family of streamlines at every footprint ra- 
dius To. See Figure 1 in Fukue & Okada (1990) for a plot 
comparing their streamlines for various values of c with the 
self-similar streamlines found by Takahara et al. (1989). By 
observing that the concavity or convexity of the streamlines 
in Figure 2 is pronounced only near the midplane, we can 
conclude that using equation (2.12) would not be an im- 
provement over our choice of geometry, as equation (2.12) 
does not capture this feature. 

Nevertheless, the findings of Fukue & Okada (1990) pro- 
vide a physical argument supporting our choice of geometry. 
Examination of equation (2.12) reveals that close to the disc 
midplane, their streamline function is indeed just a straight 
line.^ In other words, balancing the gravitational and cen- 
trifugal forces perpendicular to the flow implies straight 



^ To see this, note that to first order in z, when z « To, r ^ x 
and the right hand side of equation (2.12) is ro + cz. Hence, 
X = To + cz, which is just our s-coordinate, provided we identify 
the constant as c = 1/ tan j. Thus, c ^ 2 implies i ^ tan~^(.5) = 
26.6°. 



streamlines close to the x-axis. Equivalently, conical stream- 
lines define the path of minimum effective potential near 
the disc midplane. Moreover, Fukue & Okada showed that 
streamlines curve back on themselves (and intersect the z- 
axis) if c < 2. For streamlines to extend to infinity, it is 
required that c ^ 2. This is the requirement that the inclina- 
tion angle i < 27° (see footnote 6), which is approximately 
the opening angle of the self-similar streamlines obtained 
by Luketic et al. (2010) and shown in Figure 2. Note, how- 
ever, that the square root in equation (2.12) spoils the self- 
similarity that the streamlines would possess if the radical 
were zero. 

One would expect that a model featuring streamline 
curvature in the {x, ^;)-plaiie would lead to significantly dif- 
ferent wind solutions if the area function A{1) directly de- 
termined the critical point location. This is not the case, 
however, as the well known rocket-nozzle analogy in stellar 
wind theory revealed that (the effective) gravity, more so 
than the flow tube area, plays the role of the converging- 
diverging nozzle to facilitate a transition from subsonic to 
supersonic flow (e.g., Lamers & Cassinelli 1999; we define 
the equivalent nozzle function in §4.2). To stress this point, 
consider how the area term enters the equation of motion per 
elimination of the density gradient (by taking a logarithmic 
derivative of the continuity equation), 

1 dp _ Idv 1 dA 
p dl V dl A dl ' 

The first term on the right hand side exemplifies the out- 
come of adopting a fluid treatment: the density gradient 
(and hence the pressure gradient for a polytropic EoS) it- 
self depends on the flow acceleration. It is this term that 
gives rise to a singularity upon solving the equation of mo- 
tion for dv/dl. In turn, the second term, d\nA/dl, which 
is more a measure of streamline divergence than of the area 



(2.13) 
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between streamlines, influences the position of this singular- 
ity - the location of the critical point. This location would 
not change by much had we analytically modeled the exact 
area in Figure 2, as azimutlial streamline divergence is well 
accounted for using straight streamlines. 



3 HYDRODYNAMIC FORMULATION 

In this section wc present a general formalism for solving 
ID thermally driven wind equations under either spherical 
or axial symmetry. We adopt all of the simplifications of the 
classic Bondi and Parker problems, namely we consider the 
hydrodynamic limit, assuming a single-fluid treatment and 
inviscid, barotropic flow. Imposing these restrictions allows 
the vector Eulerian momentum equation to be integrated 
and the problem solved using a simple Bernoulli function, 
constant along a given streamline. The forces acting on a 
fluid element arc: the force of gravity from a central source, 
g£is pressure, and the centrifugal force when there is nonzero 
rotation. We use the conventional polytropic EoS in lieu of 
the energy equation. We only relax the assumption of spher- 
ical symmetry by adding an azimuthal velocity component, 
conserving the specific angular momentum of the fiuid, and 
we allow for arbitrary amounts of streamline divergence. 

Our formulation is an extension of the classic isothermal 
and polytropic Parker problems into cylindrical symmetry. 
Because of the equivalence of wind and accretion equations, 
our problem is also a generalization of the classic Bondi 
problem (Bondi 1952). Bondi's analysis entailed applying 
boundary conditions at infinity, where both the velocity and 
gravitational potential vanish. Our solution allows boundary 
conditions to be applied at any finite distance away from the 
central gravitating object. Although we do not address the 
generalized Bondi problem because our focus is on winds, it 
should be kept in mind that any explicit reference to bound- 
ary conditions taken at 'the base' - be it the disc midplane 
or the surface of the central object - can equally well de- 
note 'outer' boundary conditions appropriate to accretion 
problems. 

3.1 The Continuity Equation 

In Appendix C, we show that for the geometry of Figure 3, 
the steady-state continuity equation can be written as 

dM = p{l)v{l)A{l), (3.1) 

where dM is the differential mass-loss rate at the location 
To. The area between streamlines, A{1), can be determined 

from equation (CIO), after specifying di/dvo, the adjacent 
streamline divergence. The CIA model has no streamline 
divergence {di/dvo = 0), giving 

A{1) = 2-wdro {vo + I cos i) sin i. (3-2) 

We see by Figure 3 that the total midplane area occupied 
by the CIA model, obtained by letting / = 0, sini = 1, and 
integrating over dro from Rcia to Rout is correctly given 

as n{R'ouT ~ Rcia)- 

The Converging model, with geometry obeying tan i = 
d/vo, has streamline divergence di/dro = — cos i sin i/ro, so 

_ 2-Kdroiro + Icosif s\ni 



Notice that for the same To, both the CIA and Converging 
models have a differential base area given by Ao = A{1 = 
0) = 2-Krodro sin i. 

3.1.1 The Parker Model from the Converging Model with 
1 = 0° 

Only ratios of the area appear in the equations governing 
the flow, as in 

1 dA{l) ^ / gcosi \ 
A{1) dl \ro + lcosiJ ' ' 

where g = 1 for the CIA model and g = 2 for the Converg- 
ing model. The flow quenching factor sini does not enter 
the disc wind problem except when calculating M. We see, 
therefore, that the Converging model contains the Parker 
model as a special case, for when i = 0° , d = 0, bringing 
the converging point to the source of gravity. Then ro + 1 is 
just the spherical coordinate r, ro representing the coronal 
radius rather than the footprint distance. The only distinc- 
tion that needs to be made is that equation (3.3) formally 
does not apply in that case since A{1) = - in cylindrical 
symmetry, there is no width between streamlines because 
they all overlap on the x-axis. This can be thought of as a 
collapse to spherical symmetry, so the Parker model results, 
albeit with the adjustment that the differential base area 
becomes Ao = 2irra sin 6d6, where 6 is the spherical polar 
coordinate, instead of Ao = 27rrodro sini. 



3.1.2 Cylindrical Parker Winds: The CIA model with 
i = 0° 

Since the Parker wind solution is recovered from the Con- 
verging model at i = 0° (and with zero rotation), it is rea- 
sonable to ask if the solution to the CIA model at i = 0° 
bears any significance. It turns out that this solution was 
obtained by Skinner & Ostriker (2010) and included as a 
testbed problem in their extension of the MHD code Athena 
into cylindrical coordinates. This 'cylindrical version' of a ro- 
tating Parker wind, as they referred to it, can be viewed as a 
wind flowing perpendicular to the symmetry axis of evenly 
spaced concentric cylinders (with Ao = 2nrodz). Skinner 
& Ostriker's (2010) rotating wind test demonstrated cylin- 
drical Athena's ability to maintain steady state, transonic 
flows and conserve angular momentum in cylindrical symme- 
try. Our solutions for both the Converging and CIA models 
with i > 0° open up the possibility of allowing this test to 
incorporate the ^-dimension. 

As we demonstrate, the proper procedure for wind equa- 
tions is to take reference quantities at the footprint of a given 
streamline. The equations obtained for the CIA model at 
i = 0° by Skinner & Ostriker (2010) are seemingly the same 
Eis ours, yet the problem as they pose it is poorly formulated 
because they used reference quantities defined at infinity, 
where the pressure (and hence sound speed for a polytropic 
EoS) vanishes. Specifically, they normalized the Bernoulli 
function to c^/(7 — 1); their solutions do not suffer from 
this choice due to their assigning a value to the Bernoulli 
constant (thereby setting the location of the sonic point) a 
priori. 
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3.2 The Bernoulli Function 

By Bernoulli's theorem, the Bernoulli function is a constant 
on a streamline: 



1 2 

—V + $ + ft = constant, 



(3.5) 



for enthalpy /( = f dP/p and bulk flow velocity v — 

v'^ + v'^ + v'j = ^v{iy + v'^. We denote the Bernoulli 

constant as Bo and emphasize that, while it is to be evalu- 
ated at the boundary, 



Bo 



+ $ + h 



(3.6) 



it is a priori unknown because v{l = 0) is unknown. 

To define Bo, both the temperature T and density p 

at the base of every streamline must be specified. We take 
these quantities to be To and po, respectively. For an ideal 
gas EoS, this is equivalent to specifying the pressure at 
at every footprint location Tq. The barotropic assumption, 
dP = {dP/dp)dp, is satisfied by an ideal gas EoS combined 
with the polytropic fluid relation, T — To{p/ Po)'^~^ ■ Explic- 
itly we have P — pokTo(p/ po^ / P-nrip, which gives 



J pdp pl-^ 7 ^ 



ci \n{p/po) if 7 = 1 

c?/(7-l) if7>l, 
(3.7) 

provided that we absorb into Bo the constant term cl In(pc) 
for 7 = 1, as well as the constants of integration. 

For rotational motion in a plane under a central force, 
the specific angular momentum L is a constant of the mo- 
tion; L = xv^, where x = To + Icosi. We will present our 
disc wind results for a disc rotating at Keplerian veloci- 
ties, in which the disc angular velocity at any location ro 
is 0,K = s/GM^ro/rl. However, for treating Parker winds 
we follow Keppens & Goedbloed (1999) in allowing for ar- 
bitrary rotation rates Q, parametrizing on the equatorial 
plane by some factor of the adiabatic sound speed at the 
base, i.e. v^{l = 0) = firo = Cco, giving 



= Q.ro {—) = Cco I ^ 

\ X / \ro + I' 



(3.8) 



Keplerian rotation corresponds to (" = sJGM^i\^ = \f\,. 
Rotation therefore enters the problem as an effective poten- 
tial, 



$ $e// = — + ^ 



a2 2 / \ 2 

(, Co I To 



\/ri -1-/^-1- 2lro cos i 2 \ro -|-/ cos i 



(3.9) 



For 7 > 1 then, the Bernoulli function reads 



B„ = l^;(0' + *e//(0 + ^- (3.10) 
We treat the isothermal (7 = 1) case in §3.7. 



3.3 The Equation of Motion 

The equation of motion, which we will sometimes refer to 
as F(/, V, dv/dl) = 0, contains the velocity gradient. Critical 



points arise whenever the velocity gradient becomes unde- 
fined, i.e. when dv/dl — 0/0. Hence, the Bernoulli function 
must be accompanied by F{1, v, dv/dl) = to seek out these 
critical points. F{1, v, dv/dl) = is found by first differenti- 
ating equation (3.10) to give 



dBo dv , d$e// dp „ 
dl dl dl p dl 



(3.11) 



and then by eliminating the density gradient using the con- 
tinuity equation. The relevant derivative is given in equation 
(2.13). Further dividing by gives 

^dA 
A~dl 

(3.12) 



F{l,v, dv/dl) ; 



V dv ^ 1 d^eff 



ci dl 



dl 



0. 



3.4 Dimensionless Formulation 

Our disc wind problem depends on a total of three param- 
eters, namely, Ac, 7, and i. We find it natural to normalize 
distances to the gravitational radius, 

rg = Xoro=^^. (3.13) 

Co 

Justification for this choice is obtained by shifting one's 
viewpoint to consider all subscripts 'o' as standing for 'outer' 
rather than midplane boundary conditions. Then in the limit 
Co — >^ Coo, Tg is the so-called Bondi length. Since we will 
discuss results for both spherical winds and disc winds, we 
will differentiate disc wind bulk velocities by iiornializiiig to 
Vesc = \J GMt/fo, the escape velocity from a thin Keplerian 
disc at a distance Vo along the disc, instead of Wgac = v^Vesc- 
The HEP has several different guises in terms of these char- 
acteristic quantities, namely 



Ao = 



esc 
.2 



2 

esc 



2cl 



ro 



(3.14) 



3.4.1 Unknown Critical Point Quantities 

We introduce a quantity analogous to Ao, defined by 



* CSC 



where Cs(lc) is the sound speed at the critical point (sub- 
scripts 'c' will be used to denote quantities evaluated at 
the sonic point throughout). Since Cs(Zc) is in a one-to-one 
relationship with the critical point distance Zc, Ac is a cen- 
tral unknown. In general, Ac can only be solved for rmmeri- 
cally. Many quantities of interest such as the mass loss rate, 
the initial velocity, and the terminal velocity can be sim- 
ply expressed in terms of the ratio Ac/Ao = To/Tc. Note 
that Ac = Ao in the isothermal case. For 1 < 7 ^ 5/3, the 
ratio \c/\o is equal (by construction) to the fundamental 
constants of the problem, 

Ao ^ {Bq/cD 

Ac Cc 

where ec — ec(Xc) is the critical point constant, and is for- 
mally given by 



(3.16) 



Bo 



(3.17) 



We emphasize that Cc is a constant determined indepen- 
dently of Bo (see §3.5.1). 
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3.4-2 Dimensionless Equations 

We now rewrite the governing equations into a form suitable 
for numerical implementation. We begin by introducing the 
following dimensionless variables: 



distance along a streamline: x 



specific kinetic energy: y = 



I 



1 



2c.(Xc)2' 



Mach number: M = — , 

Cs 

Mach number squared: w = jW' 
sound speed squared: s = 



Here, Xc = Ic/vg is the dimensionless critical point distance. 
The variables y, w, and s are related by 

sw 



y = 



(3.18) 



We prefer simply keep p/po and A/Ao instead of renaming 
the density and flow tube area. Then the continuity equation 
becomes 

(3.19) 

where the mass flux density, rh, is defined as 

dM 



Ao 



(3.20) 



Note that rh = poVo — poCoMo, where Mo is the initial 
Mach number. The polytropic relation is now 



p \ _ I P^ 



Po 



Ao \Po 



(3.21) 



where by construction So = s(x = 0) = \c/\o- 

Dividing equation (3.10) by the unknown quantity 
Cs(xc)^ gives the dimensionless Bernoulli function, 



ec = y + ^Uess^ ts, 

Ao 7—1 



(3.22) 



where 

^e// = U + Ucentrif 



1 1 /_ 

Vx' + 2x/cosi + /2 ^ 2 [f 



a 



+ X cos i 



(3.23) 



and we have introduced the following quantities: 



gravitational potential: U = 



centrifugal potential: Ucentrif = 



inverse HEP: / = 



Ao' 



The above equations take their simplest form by eliminating 
any reference to po and Ao. Expressed this way, the equa- 
tions depend only on the ratio \c/\o and values taken at the 
critical point, making it clear that the accretion equations 
(with outer boundary condition Ao — >^ at vlo — > 00 but with 
finite Ac/Ao and Ac) are identical to the wind equations. 
First we square equation (3.19) and express it in terms of 



w: = sw{poCoY{\c/\o){A/ AoY{p/ Pof ■ The polytropic 
EoS, equation (3.21), permits substitution for p/po- Further 



evaluating rh at the critical point where Sc 



1 yields 



the combined continuity equation/polytropic EoS in terms 



of rh„ 



Defining 



m = m^,-T-r-ws7-i . 



A = 



equation (3.24) becomes 



,2 7+1 
A = — — WS7-1 . 
A? 



(3.24) 



(3.25) 



(3.26) 



Rewriting the dimensionless Bernoulli function in terms of 
s and w gives 



sw Ac 



7- 



(3.27) 



and equations (3.26) and (3.27) together comprise an alge- 
braic system of two equations for the two unknowns s and 
w. Once Be is evaluated, an explicit solution for w can be 
found. 

The equation of motion, equation (3.12), becomes 



1 dy 1 Ac dUcf f 1 rfA 
^ Ad^ 



0. (3.28) 



2yJ sdx s Ao dx 

This can be further simplified by letting y' = dy/dx, A' = 
dA/dx, and by defining the (minus of the) effective gravitar 
tional force as 

(C/)'cosi 



g = 



dU, 



eff 



X + fcosi 



dx 



{X^ +2x fcosi + P)i if + xcosi)^' 

(3.29) 

-F'(X) 2/) 2/') = now reads 



3.5 Critical Point Conditions 

3.5.1 The Critical Point Constant 

The value of the critical point constant Cc is found from 
equation (3.27) evaluated at the critical point. 



-2+Yo^^ff + T- 



Again since Sc = Wc = 1, we have that 



(3.31) 



(3.32) 



In general, therefore, Cc depends on the critical point dis- 
tance Xc- It is easily seen that despite the fact that the 
classic Bondi and Parker problems can have very different 
Bernoulli constants Bo, they both have the same value of 
Cc = Bo/cs{xc)^- The singular nature of their equations at 
the critical point are identical. In that spherically symmetric 
case, Cc is independent of Xc and Ac: 

1 /5-37' 



7- 



(for spherical symmetry only). (3.33) 
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Rotation breaks this equivalence because then Cc = ec(xc), 
and Xc depends on the boundary conditions. 

3.5.2 The Singularity and Regularity Conditions 

The singularity condition identifies all points at which the 
flow acceleration is undefined, i.e. all values of x for which 
F{x,y,y') = is independent of y': 

dF 

^=0. (3.34) 

From equation (3.30), the set of possible points picked out 
by equation (3.34) are those that satisfy y = s/2 (or in 
physical units, v = Cs) at Xc The regularity condition, in 
turn, defines the acceleration at this point as the slope of 
F in the (x, J/)-plane: y' = —{dF/dx)/{dF/dy), or as it is 
more commonly stated. 



dF ,dF 

ox oy 



(3.35) 



Equation 3.35 is formally derived by ensuring that 
dF{x,y,y')/dx ~ all along the solution curve, which is 
equivalent to requiring a finite jerk, i.e. that y" is bounded at 
Xc (Lamers & CassincUi 1999, §8.7). The role of the regular- 
ity condition is to ensure the continuity of the solution at the 
critical point. Since this point coincides with the sonic point 
for our problem, it marks the region where the flow loses 
communication with what is happening downstream. In the 
neighborhood of this point then, there could potentially be 
thermodynamically different situations, which would result 
in a shock - a discontinuity in y' . Physically, therefore, the 
regularity condition prevents shocks, i.e. it demands that 
nothing special happens with the flow at the critical point. 
For more complicated equations of motion, e.g. with line- 
driving included, explicit use of the regularity condition is 
required to determine the location of the critical point (e.g., 
CAK, Tsinganos et al. 1996). For thermally driven winds, it 
is not needed, but we will make use of it in §5.3 to interpret 
the negative root of the isothermal critical point equation. 

3.5.3 The Relation Between Ac and Xc 

The singularity condition combined with the equation of mo- 
tion yields a relationship between Xc and Ac; it does not 

directly determine the location of Xc (except for the isother- 
mal case when Ac = Ao). With y = s/2 at the critical point, 
equation (3.30) gives. 



Ao 

Ac 



A' 



(3.36) 



Here, Ac/Aj, = {f + XcCosi)''/(gcosi), where 5 = 1 for the 
CIA model and g = 2 for the Parker (i = 0° ) and Converging 
models. 



3.5.4 The Location of the Critical Point (s) 

The innocuous looking equation Cc = {\c/\o)Bo/(?o com- 
bined with equation (3.36) determines the location of the 
critical point. Equivalently, we can evaluate equation (3.27) 
at the lower boundary, 

Ac \Wo 



(3.37) 



where, using the definitions of y and s, we have factored 
out So = Ac/Ao. A relation between Wo and Ac follows from 
equation (3.26): 



Wo = A 



(3.38) 



With equations (3.38) and (3.32) both substituted into equa- 
tion (3.37), and noting that Ueff,o ~ — Ao + C^/2, the general 
equation that must be satisfied by a critical point is 



Ao 
Ac 

1 
2 



7 + 1 



7 

Ac 



+ ^U^ff iXc) 

'^o 



Ao + 



7-1 



(3.39) 
(3.40) 



All appearances of \o/\c in equation (3.40) are to be elimi- 
nated using equation (3.36). The resulting equation can only 
be solved numerically - with a root finder capable of detect- 
ing multiple roots - except for the classic Bondi problem. 
We solved equation (3.40) for all of its roots using simple 
bracketing and bisection (Press et al. 1992) with a tolerance 
of 10~^^ and a bracket spacing of Ax = 0.005 (although 
a spacing of 10~^ was required for 7 ^ 1.1). We observed 
that equation (3.40) almost always possess two roots. For 
the classic Parker problem, the transonic inflow solution of 
the second root satisfies equation (3.38) - there are never 
two outfiow solutions for the same set of parameters (Ao, 7). 
This issue is taken up again in §3.6.2 (and see Appendix D). 



3.6 The Polytropic Fluid Solution 

Rearranging equation (3.27) to isolate s gives 



ff- 



(3.41) 



An explicit solution to the problem, i.e. a solution comprised 
of separated functions of the dependent and independent 
variables, follows from substituting equation (3.26) solved 
for s into equation (3.41) and multiplying both sides by 

(AAc/A)-'^: 



F{w) = A-'^X(x), 



(3.42) 



where 



- 1 



'^+1 / Ac 

■ - T-Ueff 
Ao 



Note that an explicit solution cannot be found in terms of 
the kinetic energy {y = sw/2), as ridding the left hand side 
of equation (3.41) of x-dependence to give F = F{y) is not 
possible. Also note that insofar as there are solutions for 
which the ratios Ac/Ao and A/Ac are the same regardless 
of whether inner our outer boundary conditions will be ap- 
plied, the explicit solution is completely independent of the 
boundary conditions. Hence, inner or outer boundary con- 
ditions need to be applied separately to pick out the inflow 
or outflow solutions, respectively. 
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3.6.1 The Transonic Solutions 

Bondi (1952) explained why the maximum value of A = 
mlrhc occurs for the branches that pass through the criti- 
cal point. We reiterate his logic linking the maximum mass 
flux density to the transonic solutions, since it follows from 
simple mathematical considerations. For 7 > 1, both F and 
X are well-behaved functions with minimum values, so a so- 
lution ceases to exist when the right hand side of equation 
(3.42) becomes smaller than the minimum of the left hand 

side, The product A t+i X{x) is made smallest for 

some Xmin and Amax. Therefore, 

7+1 

Ama. = . (3.43) 

It is easily seen that Fmin occurs for w = 1, i.e. kj = TOc, at 
which value Fmin = .5 (7 + 1)7(7 — 1). Evidently, X^in = 
Fmin is found at Xc, which can readily be verified because 
equation (3.36) must hold at the critical point. Thus, 

Amax =Ac = l. (3.44) 

Solutions that occur for A > 1 must accordingly have 

X > Xmin, or equivalently X < Xc or x > Xc for all x, 
to satisfy equation (3.42). These are the double- valued so- 
lutions discussed below. 

In practical terms, we have just shown that to obtain the 
transonic solutions, simply set A = 1. Of course, hindsight 
into the nature of the problem led to this convenient choice 
of A, also made by Holzer & Axford (1970). 

Aside from the solution topology, our formulation of the 
general polytropic problem is now complete. In Appendix A, 
we provide formulae to compute all other variables and 
quantities of interest from the ones already given. Wc apply 
our formalism to recover the solution of the Bondi problem 
in Appendix B. 

3.6.2 The Full Solution Topology: 4 Transonic Solutions! 

As Bondi also pointed out in his 1952 paper, A acts as 
an eigenvalue in that each value of A corresponds to a 
unique set of branches, or sets of points in the (x, w)-plane 
that correspond to the two roots of the non-linear function 

R{w,x) = F{w) - A'^^X{x) = 0. He was the first to 
show that the solution possesses an X-type topology con- 
taining both single-valued branches for A < 1 and double- 
valued branches for A > 1. For A < 1, one root defines 
a subsonic branch with w < 1 and the other a supersonic 
branch in which to > 1 everywhere. This property makes 
the root-finding procedure for to straight forward, since for 
a given x, the two roots wi and W2 are always bracketed by 
< wi < 1 and 1 < W2 < 10-I-. As A ^ 1 from slightly 
below 1, these branches approach each other, bending ever 
more toward the critical point, and finally join each other 
at the single point Xc for A = 1, thus forming the transonic 
solutions. 

For the double-valued A > 1 solutions, one root defines 
a sub-critical branch with x < Xc and the other a super- 
critical branch with X > Xc always. While these branches 
cannot represent viable wind or accretion solutions by them- 
selves, they are crucial to forming complete physical solu- 
tions. That is, these branches serve the purpose of allowing 



for shock transitions to inflows or outflows (Holzer & Axford 
1970; Theuns & David 1991 and references therein). For in- 
stance, a stable wind solution consists of a transonic A = 1 

velocity profile that matches onto a subsonic, super-critical 
branch upstream of a termination shock (Velli 1994, 2001). 

It has not been emphasized in the literature that the 
solution topology of the polytropic Parker problem differs 
from the Bondi problem in one important respect: there are 
two sets of critical point solutions to choose from for every 
value of the adiabatic index in the range 1 < 7 < 5/3 (see 
Appendix D). As wc already mentioned, equation (3.40) al- 
most always possesses two roots and each one yields a family 
of solutions with the topology described above. Ultimately, 
therefore, there are two transonic outflow solutions to choose 
from for a given set of parameters. 

This second set of solutions does not arise if the bound- 
ary is at infinity, as in the Bondi problem, nor does it arise 
in the isothermal case because as 7 — >■ 1, these two sets 
of solutions approach each other and coincide for 7 = 1. In 
Appendix D, wc examine the properties of this second set of 
solutions for the classic Parker problem. Suffice it to say here 
that without rotation, only one of the two outflow solutions 
is a sought after wind solution satisfying p(x = 0)/po = 1- 

We should point out that Keppens & Goedblocd (1999) 
alluded to the presence of multiple roots of the rotating 
Parker wind critical point equation, but they incorrectly 
identified one fast rotating solution as passing through two 
critical points. To clarify our notion of two critical points, 
we reproduce their test solution in Figure 4. The top panel 
shows the two wind velocity profiles obtained for the two 
critical points; the bold outflow curve with a sonic point at 
rc ~ 3ro is the solution shown in Figure 1 of Keppens & 
Goedbloed (1999). As we explain in §5.3, what they called 
a second critical point is likely the location of the velocity 
minimum. 

In the middle panel, we show all four transonic solu- 
tions, separated into their subsonic (dashed) and supersonic 
(solid) branches. As shown in the bottom panel, only the 
bolded outflow curve satisfies the density boundary condi- 
tion p(r — >■ To) = po, and so the other must be rejected. The 
density boundary condition is also satisfied by the inflow 
curve of the second set of transonic solutions. It is tempting 
to suppose that this is always the case based on a sym- 
metry argument: owing to the mathematical equivalence of 
the wind and accretion equations, why should the density 
boundary condition preferentially be satisfied by the outflow 
solution alone? Indeed, this argument holds under spherical 
symmetry. However, one of our more interesting results is 
that the inclusion of rotation can permit both outflow so- 
lutions to satisfy the density boundary condition. How this 
comes about is addressed at length in §5.4, but we mention 
it here to draw attention to a very similar finding made by 
Cure (2004), who reported that line-driven wind equations 
also have an 'always present', second set of critical points 
when rotation is added to the problem. As in our case, he 
showed that the second family of solutions does not ordinar- 
ily satisfy the density boundary condition but can at high 
enough rotation rates. 
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Figure 4. Top panel: bulk velocity profiles of the two transonic 
solutions corresponding to the two roots of equation (3.40) for the 
parameter set chosen for a Parker wind with rotation by Kcppens 
& Gocdbloed (1999): Xo = 5.44995, 7 = 1.13, and C = 1-9. Mid- 
dle and bottom panels: Mach number and density profiles of the 
four transonic solutions corresponding to tiic two roots. The bold 
outflow solution was displayed in Figure 1 of Kcppens & Gocd- 
bloed (1999); it has rc = 2.9894 ro, Ac = 7.6807, ec = 5.9077, 
Mo = 0.5376, Vo = 0.1628 Vesc, and Voo = 0.8770i)esc. Notice 
also that it displays a Mach number minimum at r = 1.0360 r-Q. 
(The bulk velocity profile has a minimum at r = 1.0877 ro.) The 
second outflow solution has similar properties (rc = 2.3256 ro, 
Ac = 6.5036, ec = 5.7940, Mo = 0.7565, Vo = 0.2291 Dc^c, and 
foo = 0.9439 tJcsc) but docs not satisfy the density boundary con- 
dition, having instead p(ro) = 1.628 po. 



3.7 The Isothermal Solution 

For an isothermal EoS, it has been pointed out that the en- 
tire problem can be compactly solved in terms of the Lam- 
bert W function (Cranmer 2004). The many properties of 
the Lambert W function and several of its uses in physics 
can be found in Valluri et al. (2000) and references therein. 
An advantage of this solution method is that the Lambert 
W function is built into Mathematica, Maple, and MAT- 
LAB. We used it to expediently survey the parameter space 
of our isothermal disc wind solutions and to calculate the 
results presented in §4.2. 



For 7 = 1, the equation of motion is reduced to 

F^(l-Ly' + 2(^g-^)=Q. (3.45) 



Recalling that g = dUeff/dx, we obtain after integration 

\nw — w = 2[Ueff — ln(A/Ao)] + constant, (3.46) 

where we have absorbed a factor of In A~^ into the constant. 
Exponentiating and multiplying by —1 gives 

-wexp(— w) 



2 / exp[t/c//] 

-1 B 



A/Ao 



(3.47) 



Here, Fij is a constant; its value is obtained by evaluating 
equation (3.47) at the critical point, giving 



Ac 
Ao 



exp 



(3.48) 



The magnitude of Fs can be quite large (~ 10^) if the crit- 
ical point is far away, since Ac/Ao measures how much the 
flow has expanded before it becomes supersonic. 

'Operating' on both sides of equation (3.47) with the 
Lambert W function isolates —w. Hence, the solution in 
terms of the Mach number M = \/w is 



M{x) = 



\ 



AF, 



exp[^c//(x)] 
' A(x)/Ao 



(3.49) 



We have re-introduced A into equation (3.49) to distin- 
guish the transonic (A = 1) solutions from the everywhere 
sub/supersonic solutions (A < 1). Here again, setting A > 1 
yields double- valued solutions. 

The location of the critical point is obtained directly 
from equation (3.36); Xc must satisfy 



9c = 



A' 



(3.50) 



For our disc wind models, Qc is given by equation (3.29) 
and Ac/Ac = qcosi/{f + Xccosi), where q = 1 for the CIA 
and model and q = 2 for the Converging model. Due to the 
non-linear dependence on cosi, Xc must again be solved for 
numerically. For the spherically symmetric Parker model, 
meanwhile, Ac/Ac = 2/(/ + x) = 2/(r/rj) (recall that rg = 
XoTo), while the gravitational force is simply g — l/{r/rg)'^. 
We immediately recover the well known sonic point distance 
Tc = rg/2 from equation (3.50). The cylindrical Parker wind 
model (i.e. the CIA model at i = 0°) discussed in §3.1.2 also 
has g = l/(r/rg)^, while A'c/Ac = l/{r/rg), giving Vc = rg. 
We can therefore in general expect the CIA model to have 
critical points approximately twice as distant from those of 
the Converging model. 

The mass loss rate is calculated from dM = mAo by 
knowing the mass flux density, m = poCoMo. From equation 
(3.49), we find that 



Mo = V-W [-(AFb)2 exp(-2Ao + C^)]. 



(3.51) 



By the definition of the Lambert W function, we can instead 
express equation (3.51) as the nonlinear relationship 



Mo = AFs exp 



M'^ 



(3.52) 
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To an excellent approximation when A^o << 1 (valid when 
Xo » 1), therefore, the mass flux density is given by 



rh = PoCoATb exp 



-Ao + 



(3.53) 



The density distribution follows immediately from the 
continuity equation, rh = CopMA/Ag: 

p{x) _ Mo/M{x) _ ATaexp [-A„ + CV2 + >lg/2] 



Po 



A(x)/A„ 



(A(x)/Ao)X(x) 



(3.54) 

where we used equation (3.52) to obtain the second equality. 
It should be verified that p{x = 0) /po = 1 upon implemen- 
tation of these equations. 

The barometric law that we quoted in §2.1 is derived 
from equation (3.54). We consider a spherically symmetric 
(C = 0) Parker wind applied to a isothermal planetary atmo- 
sphere. Taking To to be the radius of the exobase, po is the 
density at this height. Such an atmosphere can be modeled 
as a transonic Parker wind if, in the steady state, the top 
level of the atmosphere is itself moving at speeds approach- 
ing the speed of sound at the exobase. On the other hand, 
a static atmosphere, undergoing mass loss via evaporation, 
should approximately resemble a Parker wind solution with 
m « rhc (corresponding to an everywhere subsonic solu- 
tion with a finite density at infinity). Looking to equation 
(3.49), we see that for small mass-loss rates (A << 1), the 
argument of W is small and it is valid to expand W to lead- 
ing order. ^ Thus, equation (3.54) becomes 



p = Po exp [Xo{ro/r) - Xo + Ml/2] 



(3.55) 



For a strictly static atmosphere {A4o = 0), we recover a 
barometric law, but one derived from hydrodynamic (rather 
than hydrostatic) equilibrium. 



4 RESULTS 

4.1 Estimating M using ID Models 

The scaling relationships of isothermal Parker wind solu- 
tions have been used to predict global mass loss rates from 
protoplanetary discs (Adams et al. 2004; Gorti & Hollen- 
bach 2009). For instance, Gorti & HoUcnbach (2009) com- 
bine a realistic treatment of the photevaporative heating 
process in protoplanetary discs by separately calculating the 
isothermal EUV ionization front and the FUV and X-ray 
heated neutral flow surface, thereby self-consistently deter- 
mining the base density and temperature on the flow bound- 
ary. Since their radiative transfer calculation assumed hy- 
drostatic equilibrium {M.o = 0), the hydrodynamics had 
to be included separately to calculate the mass loss rate. 
The latter is obtained by integrating the mass flux den- 
sity, rh = poVo = poCoMo, via M = 4ti rhsiniro dra- 
in other words, incorporating the hydrodynamics 'by hand' 
given Co and po requires determining a nonzero value for A4o. 
Employing isothermal Parker wind solutions offer a conve- 
nient way to accomplish this, for the HEP is obtained from 
Co and the solutions are independent of po- Given a value of 
the HEP, the sonic point distance is uniquely specified by 



equation (3.50) and the corresponding initial Mach number 
is obtained from equation (3.51). 

4-1.1 Isothermal Scaling Relationships 

Using their numerical results for po and Co, Gorti & Hol- 
lenbach (2009) obtained the mass flux density as a function 
of To by utilizing the scaling relationship for rh derived by 
Adams et al. (2004). This scaling relationship is more gener- 
ally our equation (3.53) with Keplerian rotation (^ = \/A^) 
assigned and assuming transonic outflows (A = 1): 



rh = PoCoTb exp(— Ao/2). 



(4.1) 



The quantity Fs depends on the specific streamline geome- 
try and is given by equation (3.48). The term exp (—Ao/2) 
controls the mass flux density for large HEP. The exponen- 
tial dependence results from the logarithmic enthalpy term 
and acts to suppress the wind whenever when the thermal 
energy of the gas is small compared to the escape velocity, 
e.g. when the gas is deep in the potential well of the inner 
disc region and shielded from high energy photons. 

As we mentioned in §2.3, BMS83 found the same ex- 
ponential dependence in their isothermal wind region E, 
namely rh oc exp {—Tg/2Tic), where Tg is the 'escape tem- 
perature' defined by kTg = fmipV^sc- To make the compar- 
ison explicit, we must recall equation (2.9), which says for 
7 = 1 and ^ = ro/Ric = Tic /Tg that the HEP is sim- 
ply Ao = Tg/To- For a tightly bound corona heated to the 
Compton temperature, we have that To = Tic < Tg, so we 
can indeed identify the HEP as being equal to Tg/Tic- 



4-1-2 Polytropic Scaling Relationships 

Similar agreement can be found using polytropic models, 
which are able to sample a larger range of thermodynamic 
conditions and can therefore lead to more accurate disc 
dispersal time-scale estimates. The procedure for incorpo- 
rating the hydrodynamics 'by hand' using Parker wind or 
Parker-like disc wind models is the same as that given 
above, the only change being that the initial Mach num- 

bcr is now given by Mo = (Ac/Ao)(Ao/Ac) ^(''"i) . (Corre- 
spondence with the isothermal result can be obtained using 
pc/ Po = {Xo/XcY^^'^~^^ -) Thus, for polytropic winds, the ex- 
ponential terms in equation (4.1) are replaced by a strong 
functional dependence on the temperature at the sonic point 
(Ao/Ac = Tc/To), so the mass flux density scales as 



rh = poCo{Ac/Ao) 



V -V + l 



(4.2) 



For 7 = 5/3, we recover the temperature dependence found 
by BMS83 in their Region C, namely rh oc {To/Tof - (To 
make the comparison with BMS83, Tc is to be associated 
with their 'characteristic' temperature T^h and To with Tg, 
obeying T^h < Tg < Tic-) 

By equation (3.36), we can instead express the initial 
Mach number as 



* .1 Ac ( Ac 

= a; [^^k 



T + l 
2(7-1) 



(4.3) 



W{x) : 



We see that A^o now depends sensitively on the effective 
gravitational force instead of exp {—Ucff), as well as on the 
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Figure 5. Equivalent nozzle functions N{x) for a- Keplerian 
Parker wind (topmost dashed-dotted curve), the i = 60° Con- 
verging (dashed) and CIA (solid) models, and the spherically 
symmetric Parker wind (bottom dotted curve). Rotating these 
nozzle functions about the x-a^is sweeps out the area of the de 
Laval Nozzle having steady-state flow properties identical to that 
of the wind. N{x) is normalized so that N{xc) = exp(— Ao/2), 
which is approximately Mo- The horizontal lines give the exact 
value of Mo, calculated using equation (3.51). They are, from top 
to bottom. Mo = [0.497,0.134,0.036,0.002], with corresponding 
critical points (throat positions) Xc = [0.290,0.333,0.779,0.409]. 
All nozzle functions were calculated with Ao = 11. The hump 
on the topmost curve is a bulk velocity minimum, located at 
X = 0.029. 



ratio of the flow tube area and the streamline divergence at 
the critical point. 



4.2 Isothermal Flow Properties using Nozzle 
Functions 

A useful means for gauging how the flow properties (i.e. sonic 
point distance, initial Mach number, mass flux density, and 
acceleration) of our disc wind models compare to those of 
Parker winds with or without Keplerian rotation is obtained 
by defining the equivalent nozzle function (see Parker 1963 
for the polytropic nozzle function). The well-known equsr 
tion of motion for the isothermal de Laval nozzle is, in our 
notation, 



w' = 2 



N '■ 



(4.4) 



where A'^ = A'^(x) is the cross-sectional area of the nozzle 
at any distance x- Making the comparison with equation 
(3.45), A''(x) is obtained by solving dlnN = dln^l — dUeff, 
giving 



Nix) 



Ajx) 
Ao 



exp 



-t/e//(x)- A„ 



2 2 



(4.5) 



We have normalized N{x) so that N{xc} = Ts exp(— Ao/2), 
which is the leading order approximation to M.o via an ex- 
pansion of equation (3.51). Since the critical point occurs at 
the 'throat' of the nozzle, where the cross-sectional area is 
a minimum, N{x) is a visual tool that can be used to flnd 
both the sonic point and the approximate initial Mach num- 
ber (and hence mass-flux density) by inspection. In Figure 5, 



we plot the nozzle functions of the Keplerian rotating Parker 
wind (dashed-dotted line), the Converging Model at 60° 
(long-dashed line), the CIA model at 60° (solid line), and the 
spherically symmetric Parker wind (dotted line) for Ao = 11. 
The horizontal lines correspond to the exact value of the 
initial Mach number from equation (3.51). We see that all 
of the horizontal lines except that of the Keplerian Parker 
wind with Mo ~ 0.5 intersect very near the minimums of 
the nozzle functions, showing that equation (4.1) is an ex- 
cellent approximation when A4o is small. 

A noticeable feature of the nozzle functions is the ini- 
tial hump close to the opening (near X — 0) for the mod- 
els undergoing Keplerian rotation. The spherically symmet- 
ric Parker wind nozzle is everywhere converging before the 
throat and diverging thereafter, thus ensuring that the flow 
will never decelerate. The presence of the humps indicates 
that the flow is entering a diverging nozzle (N' > 0), so that 
for initially subsonic flow, we must have w' < by equa- 
tion (4.4): the flow decelerates until reaching the top of the 
hump where N' = 0. The flow is still subsonic at this loca- 
tion, implying that the acceleration must be zero (w' = 0), 
i.e. the flow has reached its minimum velocity. The flow then 
proceeds to accelerate with the converging nozzle, traverse 
the sonic point at the throat where w = 1 and N' = (but 
w' 0), and continues to accelerate supersonically (w' > 
and w > 1) in the diverging region where N' > 0. 

Comparing the nozzle functions of the CIA and Con- 
verging models, it is clear that the former model has a sonic 
point about twice as distant as the latter, as expected. More 
distant sonic points imply smaller initial Mach numbers for 
a given HEP, and since Mo is a, direct gauge of the mass 
flux density, the total mass loss rate for a CIA wind will also 
be smaller in general. These differences all result from the 
halted expansion room of the CIA model, as will become 
clear in §4.4. Both winds experience a reduced centrifugal 
force at i = 60° , explaining why the Keplerian Parker wind 
has a signiflcantly higher initial Mach number. We can there- 
fore arrive at the result that the mass flux densities of our 
disc wind models are always bounded from below by that of 
the spherically symmetric Parker wind and above by that of 
the Keplerian Parker wind. 

More generally, plotting A''(x) allows one to easily infer 
the effects of altering the geometry of the flow or the effec- 
tive potential. For instance, the humps practically disappear 
by setting ^(x) = Ao- Recalling Figure 2, the streamlines 
found by Luketic et al. (2010) first originate from the disc 
midplane in a more vertical fashion before bending radially, 
implying that the area between streamlines indeed behaves 
as if A{x) ~ Ao for very small x- Therefore, it is likely that 
the velocity minimums would not occur in a model that cap- 
tures this feature, although it is worth noting that Luketic 
et al. (2010) observe non-monotonic radial velocity profiles 
in their fiducial run (see their Figure 5). 



4.3 Comparison with Hydrodynamical 
Simulations 

Having found an analytical wind solution for a geometry 
that approximates the simulation results in Figure 2, we now 
attempt to reproduce the shape of the sonic surface. Luketic 
et al. (2010) reported that in the region of self-similar flow, 
this shape is approximately a straight line given by Zc = ax 
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(in units of Ric), with a slope a ?s 1/4 for their fiducial 
run (see their Figure 1) and a m 1/3 for their isothermal 
run (unpublished) . Similar shapes and slopes can be seen in 
several of the figures of Woods et al. (1996) and Font et al. 
(2004). From Figure 3, we can equivalently specify a linear 
sonic surface as 



By definition, the sonic point of our analytic solution is the 
location satisfying 

- = AoXc. (4.7) 

To 

Therefore, a qualitative comparison can be made by assum- 
ing that the base of the wind has a constant HEP over the 
emitting region. Physically, this corresponds to a temper- 
ature profile satisfying To oc r^^ on the midplane. We see 
that the CIA model can accurately reproduce the sonic sur- 
face found by Luketic et al. (2010) if AoXc can be as small 
as {sini/a — cosi)~^ (which for i = 30° is about 0.88 for 
a = 1/4 and 1.58 for a — 1/3). Unfortunately, the product 
AoXc is always found to be greater than 2 for 1^7^ 5/3, 
showing that multidimensional effects cannot be ignored. 
Nevertheless, the sonic surface of the CIA model will still 
closely resemble that of the simulation results, having slope 
of about 40° for the smallest values of AoXc 

In order to make a more quantitative comparison with 
simulations, it would be necessary to know the variation of 
the HEP as a function of Tq. A detailed comparison along 
these lines will be presented in a follow up paper. 

4.4 Parameter Survey of Polytropic Transonic 
Disc Wind Solutions 

The parameter space of our disc wind models is all values 
of (7, Ao,i) that lead to transonic solutions. For a fixed 7 
and i, a parameter survey is encapsulated by a plot of the 
critical point distance as a function of HEP, which we will 
refer to as a critical point curve. Limiting our attention to 
the two intermediate angles i = 30° and i = 60° , we display 
in the top panels of Figure 6 critical point curves for select 
7 ranging from nearly isothermal (7 — 1.01) to adiabatic 
(7 = 5/3). Corresponding initial Macli numbers are given in 
the bottom panels. 

The shape of the 7 = 1.01 critical point curves can 
be qualitatively understood by considering the isothermal 
Parker model, in which rc/rg = Xc + 1/Ao ~ 1/2. For large 
HEP, the centrifugal term in the effective potential is small, 
and the Converging model will closely resemble the Parker 
model. Indeed, its critical point curve is approximately Xc = 
1/2 — 1/Ao, with the 1/Ao explaining the decreasing trend 
in X a-s Ao decreases. Meanwhile, we pointed out previously 
that the isothermal CIA model has sonic points about twice 
as distant as the Converging model, which can be seen from 
the top left panel. 

Transonic solutions at higher 7 require progressively 
smaller HEP values (e.g., higher temperatures) to make up 
for the energy lost to PdV work. Focusing on a fixed HEP 
value, e.g. Ao = 6, the critical points are shifted to ever larger 
distances as 7 increases. This occurs despite the fact that 
the gas will cool faster for larger 7 (thus lowering the sonic 



threshold to smaller bulk velocities) because the launching 
velocity becomes substantially smaller at larger 7. Transonic 
solutions with the least distant critical points have very high 
initial Mach numbers and delimit the edge of the parame- 
ter space. The dotted vertical lines that bound the HEP for 
large Xc - regardless of model and inclination angle - have 
the value 2/(7 — 1). For a given HEP in either model, the 
i = 60° solutions are launched with smaller initial Mach 
numbers than the i = 30° solutions. This is a simple con- 
sequence of the reduction in centrifugal force at larger in- 
clination angles. The extra rotational energy raises M.o for 
i = 30° above Mo for i = 60°. Conversely, the reduced cen- 
trifugal force for i — 60° permits critical points to extend to 
smaller HEP before reaching AI0 ~ 1. 

A peculiar feature of Figure 6 is the appearance of a 
'tail' on the critical point curves (shown in bold), signifying 
that at the lowest HEP for any given curve with 7 > 1.1 
(> 1.2 for the CIA model), there are two transonic solu- 
tions. As we mentioned in §3.6.2, this occurs when the out- 
flow solution corresponding to the second root of the critical 
point equation satisfies the density boundary condition. The 
tail of the critical point curve is generally in close proxim- 
ity to the disc, so these transonic solutions have a higher 
initial Mach number than the 'normal' solutions. As 7 in- 
creases, this tail grows in length, while the normal critical 
point curve shrinks. Only the tail extends to the right of 
the vertical lines at Ao = 2/(7 — 1); solutions on this side 
of the line have the property that their initial velocity ex- 
ceeds their terminal velocity. Notice that the normal critical 
point curve has an increasing slope, and one would intu- 
itively expect the critical jjoint to shift downstream as the 
temperature is decreased. Meanwhile, the tail displays the 
opposite behavior, so that more distant sonic points with 
lower initial Mach numbers have higher temperatures. We 
postpone an explanation of this behavior to §5. 

Our parameter survey reveals that disc winds possess 
solutions for 7 = 5/3, in contrast to Parker winds (with or 
without rotation). The top right panel shows that the Con- 
verging model has tail solutions to the right of the vertical 
line at Ao = 3 for a narrow range of HEP (3.5 < Ao < 4.0) 
for i = 60°. The bottom panel reveals that these solutions 
have initial Mach numbers close to unity. There are no so- 
lutions for the Converging model at i = 30° for 7 = 5/3 (at 
least not within / = lO^r^). Meanwhile, the CIA model has 
both normal and tail solutions at both angles, the normal 
ones permitting small initial Mach numbers. 

Overall, this parameter survey indicates that the CIA 
model, with its lack of streamline divergence, gives rise to 
what appears to be a more versatile wind. At small 7, for 
instance, the CIA model has transonic solutions that begin 
for Ao about half as small as the minimum HEP allowed for 
the Converging model. At larger 7, the critical points for 
the CIA model span a more appreciable range of HEP and 
altogether dominate for 7 = 5/3. 

These properties, which are solely due to geometric dif- 
ferences, are physically a manifestation of the rate of en- 
thalpy dissipation. Just considering the Bernoulli function, 
it is clear that the more rapidly that heat is liberated, the 
faster the flow must become to keep Bo constant. Converg- 
ing streamlines exhibit both lateral expansion due to stream- 
line divergence and azimuthal expansion as the wind cone 
widens, so the enthalpy can dissipate faster than it can with 
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Figure 6. Parameter survey for the CIA and Converging models for two inclination angles, i = 30° and i = 60°. Each curve has a 
constant 7; the polytropic indices between 7 = 1.01 and 7 = 5/3 are 7 = (1.1, 1.2, 4/3, 1.4, &1. 5). Initial Mach numbers in the lower 
panels correspond to the sonic point distances in the upper panels for the same Ao • The 'tails' of critical point curves with Xc decreasing 
{Mo increasing) with increasing Ao are shown in bold. Vertical dotted lines at Ao = (3, 4, 5, 6, 7) mark the value Ao = 2/(7 — 1); critical 
points to the right of this line correspond to transonic solutions with Vo > Voo ■ 



the CIA model. Conversely, the confined expansion imposed 
by the CIA streamline configuration allows the flow to re- 
tain more of its enthalpy as it expands, so that the flow can 
be launched at smaller Ao before the initial Mach number 
approaches unity. As 7 increases and the flow starts off with 
less enthalpy, a smaller compensatory reduction in HEP is 
required to launch transonic solutions compared to the Con- 
verging model. 

4.5 Disc Wind Acceleration Zones 

Each solution corresponding to the critical point locations in 
Figure 6 has an associated acceleration zone, I90, which we 
define as the distance where the flow reaches 90% of its ter- 
minal velocity, i.e. v(hu)) /Veac = 0.9\/2ec/Ac (see equation 
(A9)). A nearly isothermal wind has its temperature held 
constant to very large distances, so the acceleration zone 
extends far beyond a closer to adiabatic wind whose tem- 
perature falls off rapidly. For a given 7, the acceleration zone 
is also sensitive to the flow geometry. The CIA model has 
a greatly extended acceleration zone, again attributable to 
its geometrical conflnement curtailing adiabatic expansion. 
For example, Z90 ~ 10® rg for 7 = 1.1 compared to 10"* Vg for 
the Converging model. Apparently, the factor of two sep- 
aration in the critical point distance between these models 



translates into a two order of magnitude difference in the ex- 
tent of the acceleration zone! As 7 increases, the acceleration 
zone moves progressively closer to the equatorial plane, but 
there remains a two order of magnitude separation between 
our disc wind models. 

4.6 The CIA vs. the Converging Model 

Here we compare transonic solutions for our two streamline 
geometries. Recalling our Figure 3, the 2D disc wind model 
that we are attempting to explore with our ID solutions is 
comprised of two wind regions. The inner region hosts Con- 
verging streamlines beginning at i « 60° , that then diverge 
out to some distance along the midplane until the inclination 
angle coincides with that of the outer wind region, occupied 
by CIA streamlines at i ~ 30° . 

We limit our comparison to plotting transonic solutions 
for the same set of parameters (Ao,7) for either model with 
angles i = 30° and i = 60°. In Figure 7 we show the bulk 
velocity profiles within 21 /rg (upper panel), as well as the 
Mach rmmber and density profiles on a larger scale (lower 
panels range to x = 40), for parameters Ao = 10 and 7 = 1.1. 
The properties of these solutions are tabulated in Table 2. 
We note first off that the inclination angle only affects the 
transonic solutions in the subsonic flow regions. Secondly, 
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Table 2. Properties of the 7 = 1.1 transonic solutions plotted in Figure 7: 



Model, i 


Ao 


Xc 


Ac 


Be 


Mo 


Vo/Vesc 




Parker 


5. 


0.6496 


8.4960 


8.5000 


0.0690 


0.0218 


1.0002 


Con, 30° 


10. 


0.6327 


16.7032 


8.3825 


0.1922 


0.0608 


1.0018 


Con, 60° 


10. 


0.6270 


17.0044 


8.5058 


0.0649 


0.0205 


1.0002 


CIA, 30° 


10. 


1.7003 


18.9552 


9.4780 


0.0191 


0.0060 


1.0000 


CIA, 60° 


10. 


1.6270 


18.9680 


9.4841 


0.0110 


0.0035 


1.0000 



0.5 



0.1 



30^ 


1 




--- 60° 

><v ^^^^^ 

w-^*' ^^^^ 

.■V ^^^^ 









1 






Figure 7. Transonic disc wind solutions for the Converging model 
(unbolded) and CIA model (bolded) with 7 = 1.1 and Ac. = 10. 
Also shown isa7 = l.l, Ao = 5 Parker wind solution (dot- 
ted; the reduction in HEP is equivalent to making v^ac = Vesc). 
Top panel: bulk velocity profiles for the first llfrg. Middle Panel: 
Maeh number profiles. Bottom panel: Density profiles. Properties 
of these solutions are given in Table 2. 



there is an order of magnitude difference in the initial ve- 
locities of the CIA and Converging models. This has the 
following implication for kinematic models, which typically 
assume a velocity law based solely upon Ho , Voo , and a pa- 
rameter controlling the slope of the velocity: the values of Vo 
and the velocity gradient are quite sensitive to the type of 
wind geometry. Winds launched from i = 30° in either model 
have higher initial velocities due to the increased centrifugal 
force, in agreement with Figure 6. Despite having substan- 



tially different values of Vo, all four solutions tend to nearly 
the same terminal velocity because the Bernoulli constant 

Bo = should be unchanged for a similar set of wind 

parameters. Note that this is not obvious from Figure 7 be- 
cause the CIA model, with its extended acceleration zone, 
is yet to reach its terminal velocity at x = 40. 

Aside from the appearance of velocity minimums on the 
bulk velocity profiles, the Converging wind is well approxi- 
mated by a spherical Parker wind with an HEP half as great. 
(Since the gravitational binding energy is halved for gas in 
a Keplerian disc, a disc wind HEP twice that of a Parker 
wind preserves the ratio of the escape velocity to the ini- 
tial sound speed.) Meanwhile, the CIA model gives rise to a 
significantly slower wind. Wc can quantify the acceleration 
in the subsonic region using simple kinematics. The average 
value of the acceleration between the midplane and the sonic 
point is (a) = (v^—vD/ilc. Noting that Vc = Cs{lc), we have 
in terms of our tabulated quantities, 



(a) _ 1/Ac - (Vo/Veacf 



dir.- 2x. ' ^'-'^ 

from which we find that {a) con ~ 2. 9 {oJcia ^'^^ solu- 
tions in Figure 7. Equation (4.8) is also useful for comparing 
winds that undergo significant deceleration upon leaving the 
midplane (which occurs for large 7), in which case (o) will 
be negative. 

Although not noticeable, the bulk velocity profiles for 
the CIA model also have minimums close to the midplane. 
Gas rotating at Keplerian speeds must initially decelerate 
upon leaving the equatorial plane, independent of the flow 
parameters. This can be seen from equation (3.30); evalu- 
ated at X = for ('^ = Ao, we have that y' < provided 
q > (and q = 1 for the CIA model and 2 for the Converg- 
ing model). This situation results from the balance between 
centrifugal and gravitational forces on the disc midplane, so 
that the streamline divergence (dhiA/dl) controls whether 
or not the fiow can accelerate. The flow will decelerate so 
long as d\nA/dl > 0, which is always the case with stream- 
lines that are straight in the {x, z)-pane, for otherwise they 
would have to intersect. 

The bottom panel of Figure 7 shows that the density of 
the CIA model varies asymptotically as p oc whereas the 
Converging model has p oc This is a simple consequence 
of the continuity equation, the CIA model having an area 
term A oc / asymptotically. We see that the density at the 
critical point is smaller in the CIA model because the CIA 
wind remains subsonic out to distances nearly three times 
as large as the Converging wind; the flow has a larger dis- 
tance over which to expand. This reduced acceleration also 
accounts for why the CIA model has a larger critical point 
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Figure 8. Degenerate transonic disc wind solutions for the CIA 
model with 7 = 5/3; bolded solutions lie on the tail of the 7 = 5/3 
critical point curves in Figure 6. These tail solutions have sub- 
stantially different properties — sec Table 3. The bottom panel 
displays the behavior of the temperature within rg/2: the high- 
est velocity solutions undergo a dramatic increase in temperature 
in the region of deceleration (i.e. the gas is adiabatically com- 
pressed). 



constant and a smaller temperature at the critical point (re- 
call that Ac/Ao = Tc/To) - see Table 2. 

4.7 A Comparison of Degenerate Transonic 
Solutions 

Referring to to Figure 8, wc now examine two pairs of de- 
generate solutions for the CIA model with 7 = 5/3: one pair 
with Ao = 2.3 for i = 60° and another with Ao = 2.9 for 
i = 30° - see Table 3. EYom the Mach number profiles in 
the top panel, it is evident that each of these winds undergo 
marked deceleration before becoming sonic. Indeed, (a) is 
slightly less than zero for the tail solutions and only slightly 
greater than zero for the normal solutions. 

In the bottom panel we zoom in on the region of the 
middle panel with x ^ 0.5, in order to show the behavior of 
the temperature {T{x)/To = {p{x)/ Poy~^) just above the 
midplane. The temperature is shown in units appropriate 
for a photoionized disc heated to ~ 10* K; it is clear that 
the gas undergoes substantial adiabatic compression upon 



rising above the midplane, being heated by as much as 7QQK 
for the i = 60° tail solution. This behavior epitomizes the 
dilemma posed by degenerate solutions: if dust formation is 
to be taken into account, which of the two profiles are we to 
believe? 

The heating is not unique to 7 = 5/3 and can be ex- 
plained as follows. Since the area made available to the flow 
upon just rising above the midplane is roughly constant, the 
continuity equation implies m « pu. A decrease in velocity 
will thus be accompanied by a slight increase in density (and 
hence temperature), in general. The magnitude of this effect 
is dependent on the initial velocity. For the i — 60° solutions, 
Vo is only slightly below Co for the tail solution, so only a 
marginal increase in velocity is needed to reach the sonic 
point. The increase in temperature then acts to prevent the 
flow from immediately becoming sonic. Indeed, the i = 30° 
tail solution has a higher initial velocity than the i = 60° 
normal solution, but the latter wind has a closer sonic point 
because the former wind undergoes more adiabatic heating 
for X < 0.2. 

It is natural to suppose that one of the degenerate so- 
lutions is unstable. However, this might be difficult to un- 
cover analytically because each solution, being transonic, 
has an associated regularity condition. The regularity con- 
dition acts as a boundary condition to aid the stability anal- 
ysis and its mere presence may indicate stability (Velli 2001). 
Therefore, time dependent simulation may be a more appro- 
priate tool for assessing the stability of these solutions. 



5 DISCUSSION 

As it stands, we have solved and graphically analyzed the 
Eulerian equations for polytropic winds undergoing Keple- 
rian rotation and traversing the geometry of Figure 3. To 
fully uncover the new aspects of our solutions, we must re- 
visit the spherically symmetric Parker problem, as well as 
the 'rotating' Parker problem (i.e. a Parker wind following 
trajectories that conserve specific angular momentum). An- 
alytical considerations of the behavior of Parker winds at 
higher 7 constitute the basis of our discussion. Readers more 
interested in the observational implications of our results are 
referred to our summary, as this section is geared toward 
investigators interested in the mathematical properties of 
Parker winds. 



5.1 The Enthalpy Deficit Regime 

Depending on whether the sum of the enthalpy and the ef- 
fective potential energy terms at the boundary is positive 
or negative, wc define the flow as having an enthalpy sur- 
plus or deficit, respectively. By our dimcnsionlcss Bernoulli 
function, Cc = y + {\c/^o)Ueff + h/csiXcY', if the sum of 
the second and third terms is negative at x = 0, the kinetic 
energy y must decrease as the sum becomes less negative at 
X > to keep Cc constant. The flow can still be transonic 
because the sound speed (eventually) decreases faster than 
does the velocity. 

Since /i(x = 0)/c3(xc)^ = So/(7-l) and ?7e//,o = -\o+ 
C^/2, where So = Xc/^o, the defining condition for enthalpy 
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Table 3. Properties of the degenerate 7 = 5/3 transonic solutions for the CIA model plotted in Figure 8: 



i, (root) 


Ao 


Xc 


Ac 


ec 


Mo 


Vo/Vesc 


Vcx:> /Vesc 


30°, (tail) 


2.9 


2.6585 


9.7132 


0.9517 


0.6843 


0.4018 


0.4427 


30°, (normal) 


2.9 


15.3601 


46.4707 


0.9915 


0.1541 


0.0905 


0.2066 


60°, (tail) 


2.3 


0.8292 


3.3360 


1.1330 


0.9286 


0.6123 


0.8242 


60°, (normal) 


2.3 


1.2653 


4.5759 


1.0790 


0.6203 


0.4090 


0.6867 




Figure 9. Location of critical (sonic) points for the spherically symmetric Parker problem as a function of the HEP in the neighborhood 
of 7 = 3/2. The vertical, dotted line at Ao = 2 separates the enthalpy surplus and enthalpy deficit regimes. Critical points curves in the 
latter regime (those with Ao < 2) have a negative slope (so that higher temperatures lead to more distant sonic points) and correspond 
to decelerating transonic solutions, implying that Vo > t)oo- 



deficit flow places the following requirement on the HEP: 



For Keplerian velocities (C^ = Ao), this condition is simply 
Ao > 2/(7 — 1). We see that this regime is encountered only 
for 7 sufficiently larger than 1. From Figure 6, the Converg- 
ing model enters this regime for 7 ^ 4/3, while the CIA 
model only enters it for 7 ~ 5/3. 

As discussed by Holzer & Axford (1970), transonic 
winds require a positive Bernoulli constant, i.e. ec = 
Bo/c{xc)^ > 0. The Bernoulli function thereby permits an 
alternative definition of the enthalpy deficit regime, namely 
j/o > Ec- Equivalcntly, since Voc/vesc = \/ec/Ac (see equsr- 
tion (A9)), we must have that 

Vo > Voc- (5.2) 

In terms of the initial Mach number, this lower bound reads 
Mo > \/2ec(Ao/Ac). In the classic Parker problem, the 
transonic bulk velocity profile is a monotonically decreas- 
ing function of r for 7 > 3/2 (see e.g., Lamers & Cassinelli, 



1999). It follows that this class of solutions has Vo > Voo and 
therefore lies in the enthalpy deficit regime. 

The mathematical implication of there being a nonzero 
lower limit placed on Vo ot Mo naturally leads to the con- 
clusion reached by Parker (1960) that viable solar wind so- 
lutions satisfy 1 < 7 < 3/2. Only this class of solutions can 
have vanishingly small initial velocities, criteria that Parker 
imposed as a boundary condition. For the same reason, the 
class of 3/2 < 7 < 5/3 solutions were overlooked in followup 
treatments, e.g. that of Carovillano & King (1965). Adher- 
ing to the physical assumption that the gas is launched from 
highly subsonic speeds automatically excludes the enthalpy 
deficit regime. 

Even if plicnomcnologically motivated, discounting the 
enthalpy deficit regime prohibits insight into the full nature 
of the problem. Namely, we benefitted from the realization 
that Paxker winds undergo a regime change because the fiow 
must tap into its own kinetic energy to become transonic.'^ It 
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is informative to view the 7-dependence of the classic Parker 
problem from an enthalpy standpoint. 



5.2 The Role of Enthalpy: Spherically Symmetric 
Parker Winds 

This much has been emphasized in textbook treatments of 
the classic Parker problem: 7 = 3/2 is the one value for 
which all of the enthalpy is used up to lift the gas out of the 
potential weU, with none left over to supply kinetic energy 
(Lamers & Cassinelli 1999). The complete story is told by 
the dimensionless Bernoulli function, manipulated to read 

Ac 

Xo 



[^-Ao] +MI/2' 



(5.3) 



We can look upon equation (5.3) as either the temperature 
ratio To/Tc, or as representative of the critical point dis- 
tance, as rc/vo = Ac/2. Comparing the bracketed term in the 
denominator of equation (5.3) with the allowed HEP ranges 
(given in Table 1) allows us to infer the behavior of the crit- 
ical point as 7 increases. For transonic flows with a high 
energy input, meaning that 7 is close to 1, there is always 
a significant excess of enthalpy bcj^ond that used to combat 
gravity that can contribute to increasing the kinetic energy 
of the gas. Regardless of the HEP, Mo is small because the 
high energy input permits Vo « Co. The bracketed term is 
made smallest for large Aq, i.e., for lower coronal tempera- 
tures, so we have the intuitive notion that the smaller the 
coronal temperature, the more distant the sonic point and 
the smaller the initial Mach number. 

As 7 3/2 from below, progressively smaller HEP are 
required to compensate for the lower energy input into the 
wind as it expands. This is shown by curves on the right 
half of Figure 9. There remains a small enthalpy excess, so 
higher temperatures still lead to higher initial Mach numbers 
and less distant sonic points. We show representative initial 
Mach numbers at Tc/to = 3 and rc/r,, = 30 for 7 = 1.48 
to illustrate that Mo can be almost as small as when 7 is 
much less than 3/2, but only if the sonic point is very far 
away. (For comparison, for 7 = 1.1, rc/vo = 3 corresponds 
to Ao = 4.145 and has Mo = 0.185, while To/to = 30 cor- 
responds to Ao = 8.76 and has Mo = 1.5 x 10~^). Realistic 
solutions that traverse the sonic point within a few Vo there- 
fore require relatively high initial Mach numbers compared 
to those when 7 is small. 

For the special case 7 = 3/2, the bulk velocity is a 
constant at all radii. We show in Appendix D that for this 
particular 7, the HEP is confined to the single value Ao = 2, 
and the bulk velocity everywhere equals the initial sound 
speed, i.e. Mo = 1 and therefore Vc = Vo- The flow behav- 
ior for the class of accelerating solutions 1 < 7 < 3/2 can 
therefore be summarized as transitioning from an extreme 
surplus of enthalpy near 7 = 1 that permits Mo at low 
temperatures, to a scarcity of enthalpy near 7 — 3/2 that, 
in order to get even barely accelerating transonic solutions, 
requires high temperatures and therefore high initial Mach 
numbers and close sonic points. 

We then reach the high temperature regime (> 4x 10® /sT 



is safely in the enthalpy surplus regime for all 7 since both the 
potential and velocity vanish at the boundary, taken to be infinity. 



for solar parameters) of decelerating flow. This class of crit- 
ical point solutions with 3/2 < 7 < 5/3 was found by 
Dahlberg (1964) using a highly implicit formulation. Parker, 
responding to this finding, stated "it would be interesting to 
work out what conditions the solutions would fit to at the 
base of the corona where they start" (Parker 1965). Surpris- 
ingly, this never appears to have been done! 

Wc duly noted that the high temperature regime is char- 
acterized by an enthalpy deficit. As shown by the curves to 
the left of Ao = 2 in Figure 9, critical point distances vary 
oppositely with HEP than for 7 < 3/2; as the temperature 
decreases for a given 7, so docs the sonic point distance, 
while the initial Mach number increases. These features are 
easily explained. Loosely speaking, since the fiow is slowing 
down rather than speeding up, the gas can become sonic 
sooner only if the sound speed can quickly drop below the 
local magnitude of the bulk velocity, a scenario that is ex- 
pedited if Co is smaller (Ao larger) to begin with. More rig- 
orously, the bracketed term in equation (5.3) is negative, so 
lowering the temperature (increasing Ao) forces the initial 
Mach number to be higher, as only a high launching veloc- 
ity can compensate for the enthalpy deficit. As a secondary 
effect, the temperature is so high that lowering it somewhat 
contributes to increasing Mo- Smaller sonic point distances 
follow. We do note, however, by our representative values 
of Mo at Vc/vo = 3 and rc/fo = 30 for 7 — 1.48, that the 
initial Mach rmmber in this regime does not greatly exceed 
Mo for 7 < 3/2. Indeed, Mo can be much less than 1, but 
only for very distant sonic points. 



5.3 The Appearance of Velocity Minimums: 
Parker Winds with Rotation 

We can analytically investigate the effects of including angu- 
lar momentum by considering the 'rotating Parker problem', 
in which the streamlines are radial and spherically divergent 
but rotate rigidly with velocity = (,Co. Since stars differ- 
entially rotate, this solution is valid near the equatorial plane 
only. In a disc wind context, this solution can be viewed as 
a wind emanating from the edge of a flared, rigidly rotating 
disc; the differentially rotating, 'Kcplcrian Parker wind' so- 
lution used by Adams et al. (2004) and Gorti & HoUenbach 
(2009) corresponds to the special case C = \/\>- Results for 
the spherically symmetric Parker problem are recovered by 
setting C = 0. 

We begin by finding the relationship between fc and 
Ac, where we are calling f = r /rg = (r/ro)/\o- The effective 
gravitational force is g = [1 — {C^ /Xi)/f]/f'^ (see equation 
(3.29)), so we have by equation (3.36) that 



Ac 



2fc 



Ao l-(C/Ao)V^c' 
This equation is quadratic in fc with solution. 



To 



Ac 
4 



1± , 1-8 



(C/Ao) 



Ac/ Ao 



(5.4) 



(5.5) 



Only the positive root is satisfied by the location of the crit- 
ical point. Is there any meaning to the negative root? In 
the isothermal (Ac = Ao) case, in which the positive root 
of equation (5.5) directly yields the location of the critical 
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Figure 10. Location of critical (sonic) points vs. HEP for the Parker problem with rotation, where the (rigid) rotation rate is C = 1-0. 
This plot is to be compared with Figure 9 which has 7t = 1.5 and C = 0- Here, 7t = 1.4 and the enthalpy deficit regime sets in for HEP 
values to the left of the vertical dotted line at Ao = 2 + = 3, which divides transonic solutions with bulk velocity minimums from 
those without. Only the bolded sonic points have transonic solutions in the enthalpy deficit regime with Vo > Voo- 



point, the negative root gives the radius where the bulk ve- 
locity reaches its minimum value:® 



4 



1-4 1 



(5.6) 



This occurrence can be accounted for mathematically by re- 
calling the derivation leading up to equation (3.36). As we 
discussed in §3.5.2, the critical point must satisfy both the 
singularity and regularity conditions. The latter defines y' 
when y = 1/2. The former simply picks out all points for 
which the equation of motion does not depend on the accel- 
eration, a condition that is also satisfied if the acceleration 
is 0, i.e. at points where the velocity is a local maximum or 
minimum. In other words, while the location of the critical 
point must satisfy equation (5.4), there exist rotation rates 
when Ac = \o in which points where y' = 0, not critical 
points where y = 1/2 (and y' = 0/0 by equation (3.30)), 
also obey this equation. 

The equation governing the location of bulk velocity 
minimums in the polytropic case has the same form as equa- 
tion (5.6), but involves Smin = S(X = Xmin)- 



Xr. 



To 



4s„ 



Ac/ Ao 



(5.7) 



We are only concerned with the locations of local minimums, 
but equation (5.6) yields the location of local maximums also. 
Indeed, in the limit C — ^ 0, Tmin = is a velocity minimum for 
7 < 3/2 and a velocity maximum for 7 > 3/2. Note that purely 
decelerating flows, which are mathematically equivalent to having 
bulk velocity minimums at infinity (or maximums at Tc/to < 1), 
do not occur for sufficiently high rotation rates. 



In general, this location can only bo solved for numerically, 
as Smin is a function of both w(xmin) and A{xm.in)- How- 
ever, equation (5.6) still approximates this location when 
rrmn ^ 2ro if the flow is close to being isothermal, say 
7 < 1.05. In that circumstance, Smin « So = Ac/Ao, and 
equation (5.7) reduces to equation (5.6). 

Analyzing equation (5.7), velocity minimums will not 
arise unless r mini To > 1, implying the bound 



>^fl-^ 
2 \ 2 



(5.8) 



As it is, this inequality is not very insightful since Ac is un- 
known, but since the flow will adiabatically cool, we must 

also have Smin < So = Ac/Ao. Combining these two inequal- 
ities, we can conclude that for a given rotation rate, velocity 
minimums will not be present unless 

a2 



Ao < 2 -h C 



(5.9) 



This criteria follows more directly from the equation of mo- 
tion, as it is the condition for the flow to decelerate off the 
midplane, i.e. for y'{x = 0) < 0. Hence, inequality (5.9) is 
the statement that flows undergoing deceleration must be 
of very high temperature! It certainly contradicts our in- 
tuitive notion that higher temperatures give rise to winds 
with steeper positive velocity gradients. Evidently, the en- 
thalpy deficit regime must have Ao < 2-|-f^, and the appear- 
ance of velocity minimums is to the rotating Parker problem 
what purely decelerating fiow is to the spherically symmetric 
Parker problem. Recall that solutions in the latter problem 
with 7 > 3/2 still have monotonically increasing Mach num- 
ber profiles; never does the bulk velocity decrease faster than 
the sound speed, which would cause A4 to decrease. An effect 
of rotation is to allow this happen. Of course, the locations 
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of Mach number minimums and bulk velocity minimums do 
not coincide. It will be seen below that Mach number mini- 
mums imply accompanying bulk velocity minimums, but the 
converse is not true. 

Notice that for Keplerian rotation {(^ = Ao), the in- 
equality (5.9) is always satisfied, explaining why the bulk 
velocity profiles of our disc wind solutions possess minimums 
at all 7. It remains to explain the behavior of the tails of 
the critical point curves. For this, we must carefully study 
the enthalpy deficit regime of the rotating Parker problem. 

5.4 The Enthalpy Deficit Regime with Rotation 

The value 2 -|- is the critical HEP below which the bulk 

velocity must initially decelerate before becoming transonic. 
Each rotation rate must correspondingly lead to a unique 
polytropic index, -yt say, that determines the onset of en- 
thalpy deficit flow. By the inequality 5.1, 

Necessary (but not sufficient) criteria to reach the enthalpy 
deficit regime is for both Ao < 2 + and 7 > 7*. Sufficient 
criteria is, by definition, that Ao > 1/(7 — 1) -|- The 
inclusion of rotation lowers 74 below 3/2, which may seem 
counterintuitive but is a consequence of a tradeoff between 
the cntlialjiy and the effective potential. Consider two winds 
with the same total energy, but only one with a boost from 
the centrifugal force. The wind with rotational energy must 
be launched with a smaller temperature, and so it must in 
turn have a smaller enthalpy at the base {h ocT). However, a 
smaller temperature tends to increase the magnitude of the 
effective potential at the base (because Ueff,o = — Ao + C^/2 
and ^ Ao). Since the requirement for the enthalpy deficit 
regime can be restated as 5/(7 — 1) < \Ueff,o\, a larger 
|f7e//,o| allows for a smaller 7. 

Another effect of rotation is to smooth the transition 
from 'normal' critical point behavior, that with Tc/ro in- 
creasing with HEP, to enthalpy deficit behavior. This situa- 
tion is shown in Figure 10, which is the rotational analogue 
to Figure 9. In both figures, Ao = 2 + marks the critical 
HEP below which rc/vo begins to decrease as Ao increases. 
An obvious difference is that in Figure 10, two wind solu- 
tions arise for the same HEP to left of the vertical dotted 
line, reminiscent of our disc wind results. Only the (parts 
of) curves that are bolded are in the enthalpy deficit regime, 
having Ao > 1/(7- l)-|-C^/2. (The tail of the 7 = 1.41 curve 
is the first to meet this criteria.) As 7 is futher increased, 
the tail grows in length and the normal critical points begin 
to disappear altogether. 

The HEP bounds differ according to this behavior. If 
7 < 7t and there is only one transonic solution per HEP 
(i.e. if the curve has no tail), then 

2 + ^'<^°<;pT + y- (5-11) 

Whenever there are two critical points per HEP, the normal 
critical points have 

1 

Ao,mm(7,C) < Ao < 7 + (5.12) 

7 — i I 

The minimum HEP attained, Ao,min (7, C)i appears to be a 



priori unknown and is a function of both 7 and C,. The points 
lying on the tail, meanwhile, have 

(7,C)<Ao<2 + C'. (5.13) 

Since 74 is the polytropic index at which 1/(7 — 1) +0^/2 = 
2 + <^'^ , it marks where the upper bounds switch places. Only 
once 7 is high enough above 7* such that there is again 
only one critical point per HEP do we encounter the fully 
enthalpy deficit regime. In that case, the tail becomes the 
entire critical point curve and the HEP bound is 

;^ + ^<Ao<2 + C'. (5.14) 
Note that we recover the bounds reported in Table 1 when 

C = o. 

Another look at the transition from enthalpy surplus to 
deficit flow given in Figure 11. Rather than changing 7 for 
a fixed rotation rate, we set 7 = 1.46 and vary the rota- 
tion rate from from to C, = 1.5. In this way, the enthalpy 
deficit regime is gradually reached as 7* drops from 1.5 for 
^ = to 1.32 for C, = 1.5. Notice that the normal (enthalpy 
surplus) critical point curves, the first of which begins at 
Tc/ro ~ 1.02 for (" = 0, steadily shift upward to begin and 
end at higher sonic point distances as the rotation rate is 
increased. Slightly past C = !> however, the HEP range be- 
comes maximally confined by Ao,min(7, C) < Ao < 2 + (^'^, 
and only enthalpy deficit roots arc allowed. The latter make 
an appearance on the tail of the = 0.6 curve, the first 
curve to have Ao < 2 + = 2.36 and 7 > 7t = 1.4587. 

Figure 11 reveals that, at a fixed energy input (i.e. con- 
stant 7) near the enthalpy deficit regime, increasing the rota- 
tion rate steadily increases the minimum HEP for which crit- 
ical point solutions exist and pushes the lowest sonic point 
distance to higher values. That is to say, the centrifugal force 
permits transonic solutions to arise at lower temperatures 
overall. The fiow remains subsonic out to progressively larger 
radii because a balance must be struck between keeping the 
flow subsonic at these (still high) temperatures under low 
effective gravity and simultaneously supplying enough en- 
ergy to launch a transonic wind when there is little energy 
injected into the flow at these high 7. 

We have selected critical points at Fc/ro « 11 in 
Figure 11 to illustrate the effect that rotation has on the 
Mach number proflles. These critical point solutions have 
the properties listed in Table 4 and arc plotted in Figure 12. 
Curves with higher rotation rates always have larger initial 
Mach numbers, as expected. Bulk velocity minimums occur 
whenever the HEP is less than 2 -|- met by solutions (4) 
and higher, while we see that the Mach number profiles have 
minimums only for the three solutions in the enthalpy deficit 
regime, namely (5), (6) and (7). However, having Vo > Voo 
is not a necessary condition for there to be a Mach number 
minimum, only a sufficient one. Indeed, all critical points on 
the tail curves exhibit Mach number minimums at the rota- 
tion rates that we sampled. At even higher rotation rates, 
Mach number minimums occur for points on the normal 
critical point curves. The appearance of Mach number mini- 
mums for our disc wind solutions is therefore solely an effect 
of the high Keplerian rotation rate. 
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Figure 11. Location of critical (sonic) points vs. HEP for 7 = 1.46 for various rotation rates C = u^/co ranging from up to 1.5 for the 
Parker problem with equatorial rotation. The effect of rotation is to shift the enthalpy deficit flow regime to smaller 7, so that this regime 
(bold curves) is entered by holding 7 fixed and increasing We placed a vertical line at Ao = 2.36 to illustrate the tendency toward the 
upper bound 2 + f ^ . The numbered crosses have the properties tabulated in Table 4 and correspond to the transonic solutions plotted 
in Figure 12. 



Table 4. Properties of the 7 = 1.46 transonic solutions plotted in Figure 12: 



Label 


c 


7t 


Ao 




Tel To 


Ac 


ec 


Mo 


"^o / 'i^esc 


Voc 1 fJes 


(1) 


0. 


1.5 


2.137 


2.1739 


11.1673 


22.3347 


0.6739 


0.2348 


0.1136 


0.1737 


(2) 


0.25 


1.4923 


2.169 


2.2052 


10.9180 


21.8939 


0.6713 


0.2463 


0.1183 


0.1751 


(3) 


0.5 


1.4706 


2.268 


2.2989 


11.1108 


22.4442 


0.6639 


0.2690 


0.1263 


0.1720 


(4) 


0.75t 


1.4384 


2.435 


2.4551 


11.1.302 


22.7323 


0.6527 


0.3154 


0.1429 


0.1694 


(5) 


1.0 


1.4 


2.677 


2.6739 


10.9287 


22.6310 


0.6385 


0.3965 


0.1714 


0.1680 


(6) 


1.25 


1.3596 


3.009 


2.9552 


11.0326 


23.1550 


0.6245 


0.5196 


0.2118 


0.1642 


(7) 


1.5 


1.32 


3.485 


3.2989 


11.0584 


23.4881 


0.6119 


0.7441 


0.2819 


0.1614 



* (Ao)crit = 1/(7 - 1) + CV2 

tPor Ao = 2.435, there are two wind roots; the second has rdro = 6.2391, Mo = 0.4455, & Ac = 12.958 



6 SUMMARY & CONCLUSIONS 

We have studied generalized solutions of the classic Parker 
problem by treating spherical and cylindrical geometries in 
a unified fashion. Parker-like disc winds differ from Parker 
winds in that (i) varying the HEP for Parker winds sam- 
ples different coronal temperatures for a given Af, , whereas 
varying the HEP for disc winds samples different temper- 
atures as well as distances along the equatorial plane; and 
(ii) the purely decelerating wind regime for Parker winds 
with 7 > 3/2 is replaced by an initially decelerating wind 
that reaches a minimum speed and then proceeds to ac- 
celerate for I > To- We showed that the equivalent nozzle 
function can be used as a means to gauge whether or not 
the velocity is monotonic without actually finding the tran- 
sonic solutions. Meliani et al. (2004) also employed equiv- 



alent nozzle functions in their relativistic generalization of 
the polytropic Parker problem u.sing the Schwarzschild met- 
ric. There, winds were also found to accelerate for 7 > 3/2, 
but the velocity profiles were still monotonic. We showed 
velocity minima to be an effect of adding rotation. 

Our discussion of the spherically symmetric and rotat- 
ing Parker wind solutions showed that significant decelera- 
tion is associated with a flow regime characterized as having 
an enthalpy deficit. It is not inconceivable that this type 
of outflow can exist in an astrophysical setting, so it would 
be interesting to determine the spectral signatures of a de- 
celerating wind region. We tied the enthalpy deficit regime 
(i.e. the parameter space giving solutions with Voo < Vo) to 
the appearance of degenerate transonic wind solutions, and 
we further pointed out that the critical point behavior of 
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0.10 



2 4 6 8 10 12 14 

Figure 12. Transonic Mach number (top) and bullc velocity (bottom) profiles for the rotating Parker problem with 7 = 1.46. The 
numbered solutions have the properties listed in Table 4; each solution has a sonic point near Tc/to = 11. Solutions (5)-(7), bolded, are 
in the enthalpy deficit regime with Ao > 1/(7 — 1) + and Vo > Vao- In contrast to the spherically symmetric Parker problem, bulk 
velocity profiles are not monotonically decreasing in this regime. The bolded curves also have Mach number minimums, although this is 
due more to the increased rotation rate, i.e. enthalpy surplus solutions can also display Mach number minimums if C ^ 1- 



the second set of transonic solutions is remarkably similar 
to that reported by Cure (2004) in his study of isothermal 
line-driven stellar wind equations with rotation. Cure (2004) 
classified his new solutions as 'slow', since they obtain signifi- 
cantly smaller terminal velocities. Might those slow solutions 
be a different guise of the enthalpy deficit regime? 

Our main objective was to investigate the dynamical 
properties of two axially symmetric, thermally driven disc 
wind models, one with significant adjacent streamline di- 
vergence (the Converging model) and another with a com- 
plete lack thereof (the CIA model). We emphasize that de- 
tailed hydrodynamical simulations show that by taking into 
account the interactions of neighboring flow tubes, a self- 
similar streamline geometry emerges (recall Figure 2). We 
have neglected to mention elsewhere that CIA-like stream- 
lines have also boon found analytically from similarity so- 
lutions of idealized models for galactic superwinds. Both 
the self-similar solutions of Bardccn & Berger (1978), which 
took into account a gravitational potential, and those of Zi- 
rakashvili & Volk (2006), which did not involve gravity, are 
examples. 

Since the CIA and Converging models only differ by 
their respective amounts of streamline divergence, we can 



attribute the differences in the properties of their solutions 
as being solely due to geometric effects. Our results have im- 
plications for kinematic models that adopt a flow geometry 
similar to the Converging model for the purposes of comput- 
ing synthetic spectra to compare with observations. Namely, 
use of the Converging model will significantly overestimate 
the acceleration of the fiow if the true wind configuration 
more closely resembles the CIA model. The latter model, 
due to its smaller amount of streamline divergence, features 
a greatly extended acceleration zone, a more distant sonic 
surface, a shallower density and temperature falloff, and a 
smaller mass flux density than the Converging model for 
similar footprint conditions. Conversely, for a given mass 
flux density at the wind base, the CIA model will predict a 
higher density, implying that synthetic line profiles will ex- 
hibit stronger line absorption. Ultimately, larger error bars 
may need to be associated with the inferred mass-loss rate, 
as spherically diverging winds may tend to over-estimate it. 

Solving the time- dependent problem is likely to yield 
insights into the full domain of viable disc wind solutions. 
Especially considering that we found degenerate transonic 
solutions, uncovering the effects of time-dependence is a 
worthwhile task, one that we plan to undertake in a fu- 
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ture work. It is likely that one of the degenerate solutions 
is unstable, and this could be verified using hydrodynamical 
simulations. The alternative would be more exciting, how- 
ever, as it is conceivable that the timc-dcpcndcnt solution 
can settle upon both solutions under various circumstances. 
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In terms of the footprint value Eg = E{x = 0), the internal 
energy density is simply 

e: = j ■ 

Using equation (A2), we find that at the critical point. 



APPENDIX A: FORMULAE FOR THE 
POLYTROPIC SOLUTION 

To facilitate the usage of our solutions for inuncrical testing 
purposes, we first sketch the solution procedure. All quanti- 
ties given in Appendix A follow algebraically from the Mach 
number, A4 = ^/w, which can be determined numerically 
from the explicit solution, equation (3.42), for 1 < 7 ^ 5/3. 
Since there are multiple values of Xc satisfying equation 
(3.40), a rootfinder should screen for valid wind solutions 
by checking which roots satisfy A^o < 1 according to equar 
tion (3.38). For every critical point Xc, there is a unique 
value of Ac given by equation (3.36). 



Al The Bulk Velocity, Sound Speed, Density, &c 
Internal Energy Profiles 



We use as a characteristic velocity, Vesc = y^GMT/r^, in 
terms of which the HEP is given by Ao = Vesc/2Co. For disc 
wind solutions, simply make the substitution Vesc — > \/2Vesc, 
as the appropriate escape velocity for a Keplerian disc is 
Vesc = VGiVf./ro. 

The bulk velocity is obtained from the specific kinetic 
energy y = Ac(w/we3c)^/2, which can be found by eliminating 
s from (3.26) via y = sw/2, giving 



V 



4 



The sound speed is then simply Cs/v, 
definition, or explicitly from (3.26), 



sc = (u/Ve 



A- 



Ac_l_ 
A M 



(Al) 

c)/M by 

(A2) 



The above equations reduce to identities in the isothermal 
7=1 case in which Ac = Ao- 

The density follows straightforwardly from the poly- 
tropic EoS, 



"fA^J- 
\ A M 



(A3) 

where the second equality makes for an interesting compar- 
ison with the isothermal result, equation (3.54). Note that 
Pc/Po = (Ao/Ac)^''''''"^'. The temperature profile is readily 
obtained from equation (A2) but is most simply expressed 
in terms of the density profile as T/To = {p/ po)'^~^ ■ 

Finally, the internal energy density E is found from P = 
(7 — 1)E combined with = 'yP/p: 

(A4) 



1 



- 1) Po 

1 

27(7 - 1) 



A^J- 
A M 



2-1 
-1+1 



(A5) 



A2 The Mass Loss Rate, Initial Velocity, &: 
Terminal Velocity 

Prom equations (3.19) and (3.21) evaluated at the critical 
point, we find the critical mass flux density 



PoCo 



Ac /Xo\^^^ 



A„ V A. 



(A7) 



Recalling that Ao is a differential flow area, the total mass 
loss rate is found from M = J mAo, where the integral is 
taken over the wind region on the disc midplane (Ao = 
27rrodro sin i) or spherical boundary {Ao = 2iTro sin 9 d9). 
Since rh — poCoMo, the initial Mach number is Mo = 

(Ac/Ao)(Ao/Ac) , from which we get the initial velocity 
in escape speed units (co = Vesc/\/2Ao): 



(A8) 



Equations (A7) and (A8) both apply to 7 = 1 when casted 
in terms of the density using pc/po = {^o/ Xc)^^^~'~^^ ■ The 
terminal velocity is found by evaluating equation (3.22) at 
infinity, where both the effective potential and pressure (and 
hence s) vanish, giving 



(A9) 



APPENDIX B: THE BONDI PROBLEM 

It is instructive to apply our dimensionless formulation to 
the classic Bondi problem, in which the boundary conditions 
are evaluated at infinity. As mentioned at the beginning of 
§3.4 , rg rs = GM/c^ in that hmit. With r = r/rs, the 
potential is simply U = — 1/f (and so g = 1/r^), and the 
critical point must obey equation (3.36), giving 



1 Ac 

2 Ao' 



(Bl) 



In light of equation (Bl), the potential at the critical point 
can be written as J7c = — 2Ao/Ac. Inserting this into equation 
(3.32) gives 



■37 



7- 



(B2) 



Equations (Bl) and (B2) are the critical point conditions 
and are the same as those of the classic Parker problem. 

Determining an explicit expression for the location 
of the critical point requires knowing both the critical 
point constant, Cc, and the Bernoulli constant. Bo- We can 
then eliminate Ac/Ao from equation (Bl) by recalling that 
Xc/^o = ec/{Bo/Co). Equivalently, Ac/Ao can be found di- 
rectly be evaluating the dimensionless Bernoulli function. 
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equation (3.22), at infinity; we require j/o 
so 



Ac 

Xo 



5-37 



0, and 



(B3) 



Since Xc/Xo = c^oo/ Cs{rcY ^ equation (B3) gives a pre- 
determined relationship between the sound speed at the crit- 
ical point and the sound speed at the boundary, a situation 
unique to the Bondi problem. It is a consequence of let- 
ting Woo = and is contrary to wind problems in which 
Wo is nonzero and tied to the location of the critical point. 
With equation (B3) substituted into equation (Bl), we have 
fc = {^ — 37)74, which in physical units is the well known 
result 

GM 5 - 37 

The Bondi accretion rate is recovered from equation 
(A7), which in terms of the critical differential mass- loss 
rate reads 



Poo Coo \ 



(B5) 



With Ac = 2-Kf^r%s\n0d0, we obtain after substituting in 
equations (Bl) and (B3), 

dMc 



27rr^ sin OdOpooCa 



■37 



"(2b-i>J 



(B6) 



For the isothermal (7 = 1) case, fc = 1/2 and so 

the Bondi rate is given by equation (3.48) as Vb = 
(27rr|sin6'd6'/ylo)e3/V4- With dMc = PoCoAoTb, wc sec 
that there is no issue with calculating M despite Ao — >■ 00 
because there is no actual dependence on the area at the 
boundary. In other words, we formally have Vb = in the 
Bondi problem, as it is simply the initial Mach number of 
the flow (since dMc = Poof-oc^ooAloo, which is still valid at 
infinity because the product A^oMoc is finite). 

The isothermal Mach number profile is, from equation 
(3.49), 



M{r) 



-W 



-(Ae3/V4)2 



exp[— 1/f] 



and the corresponding density distribution is 

p{f) _ Ae'^/74 



(B7) 



(B8) 



po f'^Mif)' 

The X-type solution topology for various values of A can 
readily be explored by plotting equations (B7) and (B8). 



APPENDIX C: THE CONTINUITY EQUATION 

Here we derive the continuity equation appropriate for any 
cixially-symmetric streamline geometry consisting of straight 
streamlines in the {x, 2;)-plane. See Fukue (1990) for the ap- 
propriate generalization when streamlines possess curvature. 
In stellar wind equations with spherical symmetry, the con- 
tinuity equation is usually stated as 



M = 4'Kr'^p{r)v{r). 



(CI) 



The task is to arrive at the appropriate differential area 
function A{1) that yields the area between streamlines for 



the disc wind geometry of Figure 3. The continuity equation 
then takes the form 



dM = p{l)v{l)A{l). 



(C2) 



We begin with the steady state continuity equation in 
its coordinate-free, differential form, 

V-j = 0, (C3) 

where j = pv. Hence by the divergence theorem. 



V-jdV 



j • n dS, 



(C4) 



the area occupied by streamlines pointing along j that cross 
a surface S with surface normal h is equivalently obtained 
from 



i 



j • hdS = 0. 



(C5) 



The spherically symmetric case (j-n = j{r)) is easily treated 
by considering the fiow passing through two small solid an- 
gles dQi and dil^ at two different radii r2 > ri . The surface 
integral must vanish at both n and r2 separately, so 



j(ri)rfdni = j(r2)ridfi2. 



(C6) 



Since these radii are arbitrary, each side must equal a con- 
stant, dM, so that at any radius. 



dM = p{r)v{r)r dfi. 



(C7) 



To capture the entire flow area giving M, we recover equa- 
tion (CI) by integrating over all solid angles at constant r. 

Generalizing to the area contained between two wind 
cones at some height along the z-axis, we now work in cylin- 
drical coordinates inwhichn = 5, j(/)-i = j(Z) cos(7r/2— i) = 
j(/)sini, and dS = xdxd<f). Evaluating the surface integral 
in equation (C5) on a circular slice at some arbitrary height 
zi between two wind cones X2 > xi gives 



dM = f i{l)-hdS =j{l)smi{TTx'^)\'. (C8) 

J l=const '^1 

From Figure 3, we see that the integration limits arc from 
x\ = To + Icosi to X2 = xi + dvo + dri, where the dvo 
step sweeps out the increase in area from moving further 
out along the disc, while dri tracks the area swept out from 
streamline divergence alone. The latter can be related to 
di by projecting the arc length distance Idi onto the bold 
horizontal line: 

Idi 

dri = : — :. 

Sim 

The negative sign accounts for the decrease in the angle i 
farther out along the midplane, as we take di to be positive. 
After some algebra we arrive at 

l{di/dro) 



(C9) 



A(l) = 2Trdro(ro + Icosi) sin? 



(CIO) 



This formula for the differential area traversed by the flow 
between two straight neighboring streamlines with an ar- 
bitrary amount of streamline divergence was obtained by 
Feldmeier & Shiosman (1999) - see their equation [19].^ 

^ Note that this is only half the flow area needed to account for 
a biconical wind and hence to compute the actual mass loss rate. 
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APPENDIX D: ORIGIN OF THE 'TAILS' ON 
THE CRITICAL POINT CURVES 

In short, the tails on the critical point curves in Figure 6 are 
due to the existence of the second root to the critical point 
equation, equation (3.5.4). Insight into the nature of these 
roots can be gleaned from the spherically symmetric Parker 
problem, for which we can arrive at equation (3.5.4) by com- 
bining two separate relationships between the location of the 
critical point and the initial Mach number. The first is equiv- 
alent to equation [6] in Kcppcns & Goodblocd (1999), whoso 
quoted solutions we used to verify our numerical results. 
It follows from the singularity condition, equation (3.36), 
taken together with equation (3.38), the combined poly- 
tropic/continuity relation. Again calling f = r/vg, the ef- 
fective gravitational force is simply g — and a — f^, 
as in the Bondi problem. Subsituting \o/\c = '^/'^fc into 
equation (3.38) gives, since Ao = = Aq ^, 



Mo = 2-^ — \i 



5-3-r 

2(^-1) 



(Dl) 



A second relation between fc and Mo is found from equation 
(5.3): 



rc = 



(5 - 37)/4 



(7-1) [{^-\o)+Ml/2\ 



(D2) 



As a first indication that these equations, when combined, 
possess multiple critical points, consider the special case 7 = 
3/2. Equations (Dl) and (D2) reveal that 



Xo 



+ 1 



+ 1 



(D3) 



Substituting this back into equation (D2), we arrive at 



Mo 



4(Ao/2)^ 



(A„/2 + 1) [(Ao/2)' + 1] 



(D4) 



For Xo = 2, the only viable HEP value (see Figure 9), we 
immediately see that fc = rc/iXoTo) = 1/2, that is, ro/vo = 
1 and Mo = ^- There are thus no transonic wind solutions 
for 7 = 3/2. Values of Ao > 2 do, however, yield the locations 
of the second set of roots under examination. Notice that 
by equation (D4), these roots correspond to critical point 
solutions with Mo > ^■ 

The existence of two roots to the combined equations 
(Dl) and (D2) was known to Parker, as he mentioned them 
in his original published account using a polytropic EoS 
(Parker, 1960). Parker's analysis is insightful, using only 
Descartes' rule of signs, so we reproduce it in our notation. 
For 7 < 3/2, equations (Dl) and (D2) can be manipulated 
to read 



-3t 2(3-2t) 



xt 



^+1 
2t-i 



= 0. 



(D5) 



From the HEP bound for 7 < 3/2, we see that the first coef- 
ficient is always positive and so there are two sign changes, 
meaning there are always two critical points. ^'^ For 7 > 3/2, 



Of course, Descartes rule only applies to polynomials, so this 
analysis applies to the infinite number of 7 in the range 1 < 7 < 
3/2 that lead to integer exponents for the two fc terms. 



the exponent of the second term in equation (D5) is nega- 
tive, so the appropriate equation is now 

A^ 



Ao 



2(2t-3) 
: -v-1 



0. 



(D6) 



.7-1 / 2^ 
In this case, the factor 2(27 — 3)/(7 — 1) is never an integer 
in the range 3/2 < 7 < 5/3, so this line of analysis will not 
work. 

Dl Properties of the 'Inflow' Solutions 

The second class of transonic solutions has two interpreta- 
tions, one of them unphysical and the other physically ac- 
ceptable but very unrealistic. The latter case corresponds to 
a second outflow solution, in which the flow starts out su- 
personic and reaches a subsonic terminal velocity. The for- 
mer possibility is that of a transonic inflow solution obeying 
inner boundary conditions. This is clearly physically unac- 
ceptable because transonic flows are insensitive to conditions 
downstream of the critical point. Nevertheless, we choose to 
interpret these points as inflow solutions for the sake of clas- 
sification and comparison. With that choice, subsonic fiow 
resides at x more distant than the infiow critical point. 

To illustrate the behavior of these inflow roots and how 
they can transition to outflow roots when rotation is added 
to the problem, we survey the parameter space of both the 
spherically symmetric and Keplerian Parker wind models. 
Figure Dl is analogous to Figure 6, except that dashed crit- 
ical point curves depict inflow roots. The dashed curves in 
the top left panel are clearly of a different nature; they do not 
terminate at the vertical lines marking the values 1/(7 — 1), 
meaning that the sornc point can reside well past lO^fg, as in 
the Bondi problem. The bottom left panel shows that they 
all have approximately the same value of A^o (the terminal 
Mach number in this case). A low HEP for the inflow so- 
lutions is interpreted as establishing a large back pressure 
which can prevent the flow from becoming sonic until it is 
very close to the star. Notice that the inflow roots tend to 
the wind roots as 7 — ^ 1 and degenerate into one root in the 
strictly isothermal case in which the critical points are all lo- 
cated at fc = Xc + 1/Ao = 1/2. Finally, note that there is no 
regime change around 7 = 3/2 for the inflow roots. Again, 
for 7 = 3/2, Xc vs. HEP and Mo vs. HEP are a priori known 
and given by equations (D3) and (D4), respectively. 

The inflow roots undergo a marked change in behavior 
upon adding Keplerian rotation to the Parker problem, as 
shown in the right panels. The inflow and outflow curves be- 
come continuously connected, thereby accounting for the ap- 
pearance of a tail. The bolded portions of the outflow critical 
point curves would have stayed inflow curves had the density 
boundary condition p(x = 0)/Po = 1 remained satisfied by 
the transonic inflow solutions. However, an inflow transonic 
solution would have to traverse two sonic points in order 
to have a terminal Mach number less than unity. Since this 
is prohibited, a inflow root becomes an outflow root once 
Mo < 1. In other words, degenerate wind solutions arose 
based on a mathematical requirement, which begs a time- 
dependent solution to the problem. This occurrence may 
even be a further indication that inflow and outflow solu- 
tions are intimately coupled in such a way that the starting 
conditions of an accretion flow can lead to the subsequent 
onset of a wind (e.g., Blandford & Begelman 1998). 
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Figure Dl. Parameter survey for the Parker problem with and without Keplerian rotation. The outflow critical point curves in between 
7 = 1.01 and 7 = 1.6 in the left panels arc 7 = (1.1, 1.2, 4/3, & 1.4); the inflow roots also have a 7 = 1.5 critical point curve. The right 
panels, which have the same scaling as Figure 6, also have the intermediate 7 = (1.1, 1.2, 4/3, &c 1.4), but there are no solutions within 
10^ rg for 7 = 1.6. Neither model has solutions for 7 = 5/3. The vertical lines in the top left panel mark the values 1/(7 — 1) (recall 
Table 1), whereas they denote the values 2/(7 — 1) in top right panel. Bold portions of curves have the same meaning as in Figure 6. All 
plots in the main text display only outflow critical points. 
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